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Abstract 

We review quantum Monte Carlo methods for dealing with large shell model problems. These 
methods reduce the imaginary-time many-body evolution operator to a coherent superpo- 
sition of one-body evolutions in fluctuating one-body fields; the resultant path integral is 
evaluated stochastically We first discuss the motivation, formalism, and implementation of 
such Shell Model Monte Carlo (SMMC) methods. There then follows a sampler of results 
and insights obtained from a number of applications. These include the ground state and 
thermal properties of pf-sheW nuclei, the thermal and rotational behavior of rare-earth and 
7-soft nuclei, and the calculation of double beta-decay matrix elements. Finally, prospects 
for further progress in such calculations are discussed. 
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1 Introduction 



The description of nuclear structure began more than 60 years ago with the discovery of the 
neutron. Major milestones were the discovery of single-particle shells collective motion 
||, and their reconciliation However, a number of recent developments impose new and 
more stringent tests of our ability to describe nuclear structure. Heavy-ion induced reactions 
allow the study of nuclear behavior at extremes of temperature, angular momentum, or 
isospin j|. Increasingly precise experiments with electron ||, pion ||, kaon JF), ||, and 
nucleon |], [TO] beams probe new modes of excitation. As our understanding of supernovae 



11] and nucleosynthesis [O] is refined there is a corresponding need to know more precisely 



the relevant nuclear properties. And new pictures of nuclear structure such as the Interacting 
Boson Model |H| make implicit or explicit assumptions about the solutions of the underlying 
fermion problem that demand verification. 

The range and diversity of nuclear behavior (perhaps the greatest of any quantal many- 
body system) have naturally engendered a host of models. Short of a complete solution 
to the many- nucleon problem |TJ[, the interacting shell model is widely regarded as the 
most broadly capable description of low-energy nuclear structure, and the one most directly 
traceable to the fundamental many-body problem. Many studies have demonstrated that 
exact diagonalizations of shell-model Hamiltonians can accurately and consistently describe 
a wide range of nuclear properties, if the many-body basis is sufficiently large. Pioneering 
papers include work in the Op ||15|| , Os-ld [[UJ and O/7/2 |17| shells; more recent examples 
are [ 181, |l9l 20| . Unfortunately, the combinatorial scaling of the many-body space with 
the size of the single particle basis or the number of valence nucleons restricts such exact 
diagonalizations to either light nuclei or to heavier nuclei with only a few valence particles 
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The Shell Model Monte Carlo (SMMC) methods that have been developed in the past few 
years circumvent some of these difficulties while retaining the rigor, flexibility, and predictive 
power of traditional shell model calculations. These methods are based on a Monte Carlo 
evaluation of the path integral obtained by a Hubbard- Stratonovich (HS) transformation 



22] of the imaginary-time evolution operator. The many-body problem is thus reduced to a 



set of one-body problems in fluctuating auxiliary fields J23J. The method enforces the Pauli 
principle exactly, and the storage and computation time scale gently with the single-particle 
basis or the number of particles. Auxiliary field methods have been applied to condensed 
matter systems such as the Hubbard Model [21!}] , yielding important information about 
electron correlations and magnetic properties. 

A series of papers |26|, |27|, ^] have described the development of SMMC methods. 
The purpose of the present exposition is to summarize, in one document, the philosophy, 
formalism, numerical methods, and computational implementation of these techniques. A 
coherent presentation is particularly useful, as there is a growing body of experience in such 
calculations and there have been a few missteps along the way. We also offer a sampler 
of "first-of" calculations that illustrate the power and limitations of the SMMC approach, 
as well as some interesting physics; some of the results presented have not been published 
previously. 

Our presentation is organized as follows. In section 2 we review the rationale and defini- 
tion of the interacting shell model, some of the important characteristics of the two-body 
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interaction, and discuss some of the successes that direct diagonalization calculations have 
had in describing a broad range of nuclear properties. In section 3, we give an overview 
of SMMC methods, showing how the HS transformation can be used to express physically 
interesting observables as ratios of high-dimensional integrals. In section 4, we discuss the 
various ways in which a realistic shell model hamiltonian can be cast in a form suitable for 
the HS transform. Section 5 is devoted to a summary of the relevant aspects of MC quadra- 
ture. In section 6, we discuss various details of implementation: subtleties of the algorithm, 
numerical issues, and computational implementation. In section 7 we discuss the notori- 
ous sign problems and a practical method for their solution as well as the validation of the 
overall method. Section 8 is a sampler of various types of SMMC calculations (virtually all 
intractable by other methods), and in section 9 we offer a perspective on future work. A brief 
appendix is devoted to the relevant properties of independent-fermion quantum mechanics. 
Our review covers work through December, 1995. 

2 Review of the shell model 

In this section, we present a brief definition and overview of the nuclear shell model. Our 
goals are to establish conventions and notation and to give the non-expert some appreciation 
for the special features of the nuclear problem relative to other quantal many-body systems. 
More detailed discussions can be found in several texts || [ID], [33| . 

The notion of independent particles moving in a common one-body potential is central to 
our description of atoms, metals, and hadrons. It is also realized in nuclei, and the shell 
structure associated with the magic numbers was first put on a firm basis in 1949, when the 
magic numbers were explained by an harmonic oscillator spectrum with a strong, inverted 
(with respect to the atomic case) spin-orbit potential 

But nuclei differ from the other quantal systems cited above in that the residual interac- 
tion between the valence fermions is strong and so severely perturbs the naive single-particle 
picture. This interaction mixes together many different configurations to produce the true 
eigenstates and, because of its coherence, there emerge phenomena such as pairing, mod- 
ification of sum rules, deformation, and collective rotations and vibrations. An accurate 
treatment of the residual interaction is therefore essential to properly describe nuclei. 

The nuclear shell model is defined by a set of spin-orbit coupled single-particle states with 
quantum numbers Ijm denoting the orbital angular momentum (I) and the total angular 
momenta (j) and its ^-component, m. Although non-spherical one-body potentials are a 
common efficiency used in describing deformed nuclei, for the rotationally invariant hamil- 
tonians used in SMMC so far, these states have energies ey that are independent of m. The 
single-particle states and energies may be different for neutrons and protons, in which case 
it is convenient to include also the isospin component t% = ±1/2 in the state description. 
We will use the label a for the set of quantum numbers Ijm or Ijmt^, as appropriate. 

The particular single-particle states included in a given calculation depend upon the physics 
being addressed, but at least one major shell is believed to be necessary to adequately 
describe low-lying states of a given nucleus. We will use N s to denote the number of such 
states. Thus, for example, N s is (12, 20, 32, 44) for either neutrons or protons in the (IsOd, 
lpOf, 2sldg 7 / 2 h u/ 2, 2plfh 9/2 i 13 / 2 ) shells. 
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The totality of Pauli-allowed configurations of the valence nucleons in the single particle 
states defines the model space in which the many-body Hamiltonian acts. Computing the 
dimension of this space (number of such many-body states) is a straightforward combinatorial 
exercise. As noted in the Introduction, the dimension increases strongly with either N s or the 
number of valence nucleons, and can vastly exceed 10 8 for realistic applications of current 
interest. The size of the Hamiltonian matrix to be considered can be reduced somewhat 
by exploiting rotational and isospin invariance properties. Even so, constructing such a 
matrix and finding its lowest eigenvalues and eigenstates is difficult in some cases of interest 
and impossible in most, and thermal properties are completely inaccessible without ad hoc 
assumptions. 

The shell model Hamiltonian can be written in the form H = Hi + H 2 where 

Hi = ^2e a ala a , (2.1a) 

a 

H 2 = 2 E VafrsaWpasa-y ■ ( 2 - lb ) 

Here, and a are fermion creation and annihilation operators, and the V are the uncoupled 
matrix elements of the two-body interaction. These latter must respect rotational and time- 
reversal invariance, and parity conservation. To make explicit the rotational invariance and 
shell structure, we can rewrite the two-body Hamiltonian as 

* abed J M 

= J E E + ^(l + <M] V2 Vj A (ab, cd) ]T Al M {ab)A JM {cd) , (2.2) 

abed J M 

where the sum is taken over all proton and neutron single-particle orbits (denoted by a, b, c, d) 
and the pair creation and annihilation operators are given by 

Al M (ab) = Y, (jam a j b m b \JM)a) bm a) ama = -[a) a x a] b ) JM , (2.3a ) 

m a m b 

A JM (ab) = Yl Uam a j b m b \JM)a jama a jbmb = [a ja x a jb ] JM . (2.3b ) 

m a m b 



The Vj(ab, cd) are the angular-momentum coupled two-body matrix elements (TBME) of 
a scalar potential V(ri,r^) defined as 



Vj(ab,cd) = (fe(ri) x ^ jb (?2)) JM \V(n,r 2 )\[Tpdn) x ^M)\ JM ) , 



JM\ 



(2.4) 



(independent of M) while the anti-symmetrized two-body matrix elements Vj i (ab, cd) are 
given by 

Vj A {ab, cd) = [(1 + 6 ab )(l + 6 cd )]~ 1/2 [Vj(ab, cd) - (-iy^- J Vj(ab, dc)] . (2.5) 

There are a number of methods for deriving the residual interaction V (which acts in the 
model space) from the underlying "bare" internucleon interaction [J33, |35|, RB|]. Although a 
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complete specification of V requires many two-body matrix elements p7| , [38] (e.g., 63 in the 
complete s<i-shell and 195 in the complete p/-shell), successful interactions that reproduce a 
large body of experimental data show a few simple features (e.g., isoscalar pairing, attractive 
quadrupole, repulsive dipole, ...), with the rest of the TBMEs being small and random 
A venerable approximation is to truncate the interaction to just the isoscalar 



pairing and the quadrupole interaction 40 



Various aspects of the eigenstates of H can be understood through such approximations 
as Hartree-Fock [fil| , [12| or Random Phase ||43|| . However, such approaches are unsuitable 



for precise work. The preferred method is to expand the wavefunction in a many-body basis 
(truncated by some criteria, if necessary) and then construct and diagonalize the resulting 
hamiltonian matrix. Several different approaches to these tasks have been followed over 
the years. In the jj-coupling scheme, developed by the Rochester-Oak Ridge group |I7[, the 
antisymmetric iV-particle states in each j-shell are first built recursively, before the multishell 
model space is constructed from these single-shell states. The most time-consuming part of 
this method is the calculation of the single-particle coefficients of fractional parentage to 
form antisymmetric iV-particle states and it has been found that this scheme has severe 
limitations when the shell dimensions are large. In the m-scheme approach |4"4" , |4"5| each 
iV-particle Slater determinant is represented by a computer bit using the binaries and 1 to 
represent unoccupied and occupied levels. An advantage is that the m-scheme reduces the 
calculation of matrix elements (in second quantization) to very efficient logical operations. 
However, it becomes impractical when the dimensions involved get very large (i.e., with 
increasing numbers of particles and/or orbitals). The OXBASH code combines jj-coupling 
and m-scheme |46|, [47j by constructing the basis states in the m-scheme approach, but 
following the jj-coupling scheme in the diagonalization of the Hamilton matrix. Recently, 
the Drexel group has developed an alternative approach based on a pair-truncated fermion 
model space and permutation group ideas p0| . A more detailed description of the various 



shell model concepts and codes can be found in [48 



In conventional shell model applications the dimension of the model space makes a complete 
diagonalization of the Hamiltonian impractical. As one is usually interested only in the 
nuclear spectrum at low energies, the diagonalization is often performed using the Lanczos 
algorithm, which is an efficient way to find the few lowest (or highest) eigenvalues and 
eigenvectors of a large matrix [f4"9"|, pCfl. 



Physical applications of the shell model were pioneered by Lawson and by Kurath and 
Cohen in the p-shell 



30, 15 



As the dimensions of the model spaces in the s<i-shell are still 
not too large, the conventional shell model approaches are well suited and have been quite 
extensively and successfully applied. In particular, Wildenthal has performed systematic 
studies of all srf-shell nuclei and, using a specially designed effective interaction ||37|| , has 

However, 
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been able to reproduce energies and transition moments in these nuclei 
Figure 2.1 illustrates that accurate results can be obtained only in an untruncated basis. 

Beyond the sd-shell, applications of the conventional shell model become more problem- 
atic as the large dimensions of the model space make diagonalization impossible. To date, 
systematic studies of the A = 48 isobars, involving model spaces with about 10 6 Slater deter- 
minants, have been the limit for the conventional shell model approaches |51 |. Nevertheless 



the successful description of these nuclei underlines again that the shell model concept is also 
the method of choice in heavier nuclei - if such applications were only possible. As noted by 
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order 1 2 3 4 12 Exp 

dimension 1 13 261 2345 11398 93710 

FIG 2.1 Convergence of the spectrum of 28 Si with increasing dimension of the model space; 
"order 12" corresponds to the complete sci-shell space PT . 



Caurier et al. [[Hj] , it took two generations of hardware and software development to extend 
shell model calculations from A = 44 to A = 48. Thus, conventional shell model calculations 
for nuclei heavier than A = 50 appear out of reach in the near future, even assuming an 
optimistic increase of computer capabilities. 
Facing the impracticality of complete shell model calculations, Poves, Zuker, and collab- 



orators |52| have used a truncation method for the shell model in the lpOf space. When 



applied to observables like the ground state energy or Gamow- Teller strength, the conver- 
gence appears fast enough to obtain reasonable approximations to the (complete) shell model 
results after a few truncation scheme iterations. As an example, Fig. 2.2 shows the Gamow- 
Teller strengths in 54 Fe and 56 Fe calculated at various levels of truncation; there is clear 
convergence to the complete result obtained with the SMMC methods (see Section 8). 



3 Overview of Monte Carlo Methods 

In this section, we give an overview of SMMC methods. In particular, we describe what the 
methods are (and are not) capable of calculating. We then discuss the essence of the HS 
transformation, which is at the heart of the calculations. 

The SMMC methods rely on an ability to deal with the imaginary-time many-body evo- 
lution operator, exp(—(3H), where /3 is a real c-number. While this does not result in a 
complete solution to the many-body problem in the sense of giving all eigenvalues and eigen- 
states of H, it can result in much useful information. For example, the expectation value 
of some observable Cl in the grand canonical ensemble can be obtained by adding to H a 
term —fi n N — fj, p Z, (// n and \i p are the neutron and proton chemical potentials) and then 
calculating 
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FIG 2.2 Comparison of the total Gamow- Teller strengths B(GT + ) for 54 > 56 Fe in a series of 
direct diagonalizations with decreasing level of truncation ]52] with the complete results 
obtained with the SMMC method (solid symbols at t = 10) |55fl . The open symbols at 



t = 10 show the extrapolated no-truncation results of p3] . 



Here, f3 = T 1 is interpreted as the inverse of the temperature T, and the many-body trace 
is defined as 

TrX = J2(i\X\i) , (3.2) 

i 

where the sum is over all many-body states of the system. Similarly, if Pa = 5 (A — N) is the 
projector onto states with A nucleons (actually the product of separate neutron and proton 
projectors), the canonical ensemble is defined by 

Tt a X = ^\PaX\z) , (3.3) 

i 

and the associated expectation value is 

(O). = Tr ^^ . (3.4) 

In the limit of low temperature (T — >• or (3 — > oo), both the grand-canonical and canonical 
traces reduce to ground state expectation values. Alternatively, if |$) is a many-body trial 
state not orthogonal to the exact ground state, |^), then e~^ H can be used as a filter to refine 
|$) to as f3 becomes large. An observable can be calculated in this "zero temperature" 
method as 

($|e-g*fi e -4*|$) ^ (<P|fi|tt) f 5) 

($|e-/» A |$) (^|^) ' 
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If Q is the hamiltonian, then (3.5) at (3 = is the variational estimate of the energy, 
and improves as (3 increases. Of course, the efficiency of the refinement for any observable 
depends upon the degree to which |<3>) approximates \$>). 

Beyond such static properties, e~^ H allows us to obtain some information about the dy- 
namical response of the system. For operators W and Cl, the response function i?n(r) in the 
canonical ensemble is defined as 

Mr) - - (rt(T)Sl(0)) A , (3.6) 

where ftt(r) = e rH ^e' rH is the imaginary-time Heisenberg operator. As we shall see in 
section 8 below, interesting choices for Cl are the dj for particular orbitals, the Gamow- 
Teller, Ml, or quadrupole moment, etc. Inserting complete sets of A-body eigenstates of H 
|/)} with energies E i: f) shows that 

Rn(r) = |£e-^|(/|Q|i)|V^/-^, (3.7) 

if 

where Z = J2i e~^ Ei is the partition function. Thus, Rq.{j) is the Laplace transform of the 
strength function Sn(E): 

/oo 

e~ TE S n {E)dE ; (3.8a ) 
-oo 

S n (E) = We-^ifMi^HiE-Ef + E,). (3.8b) 

fi 

Hence, if we can calculate Rq(t), Sq(E) can be determined. Short of a full inversion of the 
Laplace transform (which is often numerically difficult; see section 6.4 below), the behavior of 
Rn(r) for small r gives information about the energy-weighted moments of Sq. In particular, 

S n (E)dE = -J2e-^\(f\n\i)\ 2 = (^n) A (3.9) 

-oo Z • 

is the total strength, 

/°° I 
S n (E)EdE = -^Y.e~ PE *\(f\m 2 (E f - E t ) (3.10) 
-oo Zj if 

is the first moment, 

/oo I 
S n (E)E 2 dE=-Y,e~^\(f\m 2 (E f -E i r (3.11) 
-oo Zj if 

is the second moment, and so on. (In these expressions, the prime denotes differentiation 
with respect to r.) 

It is important to note that we cannot usually obtain detailed spectroscopic information 
from SMMC calculations. Rather, we can calculate expectation values of operators in the 
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thermodynamic ensembles or the ground state. Occasionally, these can indirectly furnish 
properties of excited states. For example, if there is a collective 2 + state absorbing most 
of the E2 strength, then the centroid of the quadrupole response function will be a good 
estimate of its energy. But, in general, we are without the numerous specific excitation 
energies and wavefunctions that characterize a direct diagonalization. This is both a blessing 
and a curse. The former is that for the very large model spaces of interest, there is no way in 
which we can deal explicitly with all of the wavefunctions and excitation energies. Indeed, 
we often don't need to, as experiments only measure average nuclear properties at a given 
excitation energy. The curse is that comparison with detailed properties of specific levels is 
difficult. In this sense, the SMMC method is complementary to direct diagonalization for 
modest model spaces, but is the only method for treating very large problems. 

It remains, of course, to describe the Hubbard-Stratanovich "trick" by which e~@ H is 
managed. In broad terms, the difficult many-body evolution is replaced by a superposition 
of an infinity of tractable one-body evolutions, each in a different external field, a. Integration 
over the external fields thus reduces the many-body problem to quadrature. The following 
schematic discussion serves to illustrate the central idea; important details will be added in 
the following sections. 

With some rearrangement, as discussed in Section 4, the many-body Hamiltonian (2.1) 
can be written schematically as 

H = e6 + \v66, (3.12) 

where O is a density operator of the form a^a, V is the strength of the two-body interaction, 
and e a single-particle energy. In the full problem, there are many such quantities with 
various orbital indices that are summed over, but we omit them here for the sake of clarity. 

All of the difficulty arises from the two-body interaction, that term in H quadratic in O. 
If H were solely linear in O, we would have a one-body quantum system, which is readily 
dealt with (see the Appendix). To linearize the evolution, we employ the Gaussian identity 



-pit = J £111 [°° dae-^ 2 e-P h ; h = eO + sVaO . (3.13) 
V 2n J-oo 



Here, h is a one-body operator associated with a c-number field a, and the many-body 
evolution is obtained by integrating the one-body evolution JJ a = e~ l3h over all a with a 
Gaussian weight. The phase, s, is 1 if V < or i if V > 0. Equation (3.13) is easily verified 
by completing the square in the exponent of the integrand and then doing the integral; since 
there is only a single operator O, there is no need to worry about non-commutation. 

With an expression of the form (3.13), it is now straightforward to write observables as 
the ratio of two integrals. For example, the canonical expectation value (3.4) becomes 

{S1)A ~ Jdn-iW^u. • <3 - 14) 
which can be more conveniently written as 

(n) A =-^-; (3.15a) 
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W a = G.TtaU,; G a = e-ll^ 2 , il a = ^E^L (3.15b) 

Thus, the many-body observable is the weighted average (weight W a ) of the observable 
Q a calculated in a canonical ensemble involving only the one-body evolution U a . Similar 
expressions involving two a fields (one each for e~ rH and e~^~ T ^ H ) can be written down for 
the response function (3.6) and all are readily transcribed to the grand-canonical ensemble 
(3.1) or zero-temperature case (3.5). 

An expression of the form (3.15) has a number of attractive features. First, the problem 
has been reduced to quadrature — we need only calculate the ratio of two integrals. Second, 
all of the quantum mechanics (which appears in Q a ) is of the one-body variety, which is 
simply handled by the algebra of N s x N s matrices. The price to pay is that we must treat 
the one-body problem for all possible a fields. 

For a realistic hamiltonian, there will be many non-commuting density operators O a 
present, but we will always be able to reduce the two-body term to diagonal form. Thus for 
a general two-body interaction in a general time-reversal invariant form we write 

H = ]T (e* a 6 a + e a 6 a ) +^£K {O a , O a } , (3.16) 

a a 

where O a is the time reverse of O a . Since in general, [O a , Op] 7^ we must split the interval 
P into N t "time slices" of length A/3 = (3/N t , 

e -0H = [ e -A^]* (3 17) 

and for each time slice n = 1, . . . , N t perform a linearization similar to (3.16) using auxiliary 
fields a an : 

e -Af3H _ joo jj^ ^ d<T an dal n ^\V a \ ^ e -Af3j2jV a \\<r an \ 2 e -Aph n . 

K =E a (e a + s a V a a an )d a + (e a + s a V a a* an )6 a . (3.18) 

Note that because the various O a need not commute, (3.18) is accurate only through order 
A/3 and that the representation of e~ A/3h must be accurate through order A/5 2 to achieve 
that accuracy. 

Upon composing (3.18) many times according to (3.17), we can write expressions for 
observables as the ratio of two field integrals. For example, (3.15) becomes 

{n)A = fVaW„ ' (3 ' 19a } 

where 

W a = G a Tr A U ■ G u — e -^E a J^\M 2 . 

Q a = ; Va = [[ [[d(T an d(T an — , (3.19b ) 

±*aU n=l a \ Z7F / 

and 
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U = U Nt ...U 2 U 1 ; U n = e - A ^; 
h n = J2 ( £ a + s a V a a an ) O a + (e a + s a V a a* an ) 6 a . (3.19c ) 

a 

This is, of course, a discrete version of a path integral over a . Because there is a field variable 
for each operator at each time slice, the dimension of the integrals T>a can be very large, 
often exceeding 10 5 . Note that because of the errors in (3.18), the errors in Eqs. (3.19) are 
of order A/?, so that high accuracy requires large N t and perhaps extrapolation to N t = oo 
(A/3 = 0). 

Two steps are then necessary to implement a SMMC calculation. First, we must learn 
how to write expressions like (3.16, 3.19) for a realistic hamiltonian. Second, we must learn 
how to efficiently compute the resulting ratio of high-dimensional integrals. These tasks are, 
respectively, the subjects of the following two sections. 

4 Decomposition of the Hamiltonian 

To realize the HS transformation, the two-body parts of H must be cast as a quadratic form 



in one-body density operators O a . As discussed in Ref. [E7|, there is considerable freedom 



in doing so. In the simplest example, let us consider an individual interaction term, 

H = ala^a^as , (4-1) 

where a], are anti-commuting fermion creation and annihilation operators. 
There are then two ways to proceed: we can group (1, 3) and (2, 4) to get 

H = a\a 3 a^a^ — a{a 4 5 23 (4.2a) 

1 1 1 

= -a\a 4 5 2 ?, + -[a\a 3 , a\a^\ + -(a[a 3 + ala^) 2 - -(a\a 3 - a\a A ) 2 , (4.2b) 

or group (1,4) and (2,3) to get 



H = — a\di a\a 3 +a\a 3 524, (4.3a 
11 ■ 1 

= ^1^3^24 — 7T[ala4, «2 a 3] — T( a l a 4 + a 2 a 3) 2 + T( a l a 4 — a 2 a 3) 2 • (4.3b 



2 L 1 1 * 1 4 V " * 4 



The commutator terms in both (4.2b, 4.3b) are one-body operators, but now the quadratic 
forms are squares of density operators. We refer to Eq. (4.2b) as the 'direct' decomposition 
and Eq. (4.3b) as the 'exchange' decomposition. At the formal level, it is also possible to 
consider a pairing decomposition p7 |, although no practical calculations have been done 



with that scheme. 

For any general two-body Hamiltonian, there is the freedom to choose between the direct 
and exchange formulations and it is particularly convenient to use quadratic forms of density 
operators that respect rotational invariance, isospin symmetry, parity conservation, and the 
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shell structure of the system. Although the exact path integral result is independent of the 
scheme used, different schemes will lead to different results under certain approximations 
(e.g., mean field or static path approximation [5B|). The choice of decomposition will also 
affect the rate of convergence of numerical results as Nf — > oo, as well as the statistical 
precision of the Monte Carlo evaluation. 

We note that the two-body Hamiltonian for fermion systems is completely specified by 
the set of anti-symmetrized two-body matrix elements Vf~(ab, cd) of Eq. (2.5) that are the 
input to many standard shell model codes such as OXBASH [fi6]| . Indeed, we can add to the 
Vj~(ab, cd) any set of (unphysical) symmetric two-body matrix elements Vf(ab, cd) satisfying 

Vj s (ab, cd) = (~l) jc+jd - J Vj S (ab, dc) , (4.4) 

without altering the action of H2 on any many-fermion wave function. However, note that 
although the Vf(ab, cd) do not alter the eigenstates and eigenvalues of the full Hamiltonian, 
they can (and do) affect the character of the decomposition of H2, as is shown below. In 
what follows, we define the set of two-body matrix elements Vf(ab,cd) that may have no 
definite symmetries as 

Vj N (ab, cd) = Vj A (ab, cd) + Vf{ab, cd) , (4.5) 
allowing us to write the two-body Hamiltonian (2.2) as 

H 2 = )EE[(1 + S ab ){l + 5 cd )] 1/2 Vf(ab, cd) A\ M {ab)A JM (cd) . (4.6) 

abed J M 

To decompose H2, we perform a Pandya transformation to recouple (a, c) and (b, d) into 
density operators with definite multipolarity, 

pKAi{ab) = {jam a j b m b \KM) 

a j a m a a jbm b > (4.7) 

m a ,m b 

where dj artla = (— l) Ja+ma a Ja _ ma . Then H 2 can be rewritten as 

H 2 = H' 2 + H[- (4.8) 

H'2 = \ EE^K hd ) E(- 1 )" /J r pK-Ai(ac)p KM (bd) , (4.9) 

* abed K M 

where the particle-hole matrix elements of the interaction are 
E K (ac, bd) = (-iy*+ic Ej(-l) J (2J + 1 



ja jb J 1 
3d je K J 

Vj N (ab,cd)^(l + 5 ab )(l + 5 cd ) , (4.10) 



and H[ is a one-body operator given by 



H[ = J2z'adPoo(a,d), (4.11) 

ad 
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with 



e' ad = ~\ E E(-l) J+ja+j VJ + l)^j=fj V J N ( ab > hd )^ 1 + W + U • (4.12) 

Note that adding symmetric matrix elements is equivalent to using the exchange decompo- 
sition for some parts of the interaction. The freedom in choosing the combinations of direct 
and exchange decomposition is then embodied in the arbitrary symmetric part of the matrix 
elements . 

Introducing the shorthand notation i = (ac),j = (bd), we can write Eq. (4.9) as 

H'2 = \Y.Y. E k^3)(-^) M Pkm^)pk-mU) ■ (4.13) 

Z ij K 

Upon diagonalizing the matrix Ex(i, j) to obtain eigenvalues \x a and associated eigenvectors 
vkol-, we can represent H' 2 as 

H' 2 = W X K(a)(-l) M pKM(a)p K -M(a) , (4.14) 



2 Ko 



where 



Pkm{o) = J2PKM(i)v Ka (i) ■ (4.15) 

Finally, if we define 

Q K M{a)= ^ 2(1 ^ Mp) (pK M {a) + (-1) M Pk-m{o)) , (4.16a) 

pKM ^ a) = ~ y/2(i+s M0) (P KM ^ ~ (-1) M PK- M (a)) , (4.16b ) 



then H' 2 becomes 

H' 2 E (<&*(«) + P 2 km^)) ■ (4-17) 

Z Ka M>0 

This completes the representation of the two-body interaction as a diagonal quadratic form 
in density operators. We then couple auxiliary fields <Jkm(&) to Qkm and r KM (a) to Pkm 
in the HS transformation. (The latter are not to be confused with the "imaginary time" r.) 

In the treatment thus far, protons and neutrons were not distinguished from each other. 
Although the original Hamiltonian H 2 conserves proton and neutron numbers, we ultimately 
might deal with one-body operators PKM{(ip-,b n ) and PKui&n^p) (n,p subscripts denoting 
neutron and proton) that individually do not do so. The one-body hamiltonian h a appearing 
in the HS transformation then mixes neutrons and protons. The single-particle wavefunctions 
in a Slater determinant then contain both neutron and proton components and neutron 
and proton numbers are not conserved separately in each Monte Carlo sample; rather the 
conservation is enforced only statistically. 
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It is, of course, possible to recouple so that only density operators separately conserving 
neutron and proton numbers (pKM(a P , b p ) and pKuifln-, b n )) are present. To do so, we write 
the two-body Hamiltonian in a manifestly isospin-invariant form, 

H2 = \ E E K 1 + W + U] 1/2 Vj N T (ab, cd) E A^ MTz (ab)A JT ;MTz(cd) , (4.18) 

4 abed JT MTz 

where, similar to the previous definition (4.6), the pair operator is 

1 1 

Aj T ;MTz( ab ) = E Uam a j b m b \JM)(-t a -t b \TT z )a] bmbtb a] amata . (4.19) 

m a ,m b 

Here (|, t a ), etc. are the isospin indices with t a = — | for proton states and t a = \ for neutron 
states, and (TT Z ) are the coupled isospin quantum numbers. The two-body Hamiltonian can 
now be written solely in terms of density operators that conserve the proton and neutron 
numbers. Namely, 

H 2 = H[ + H' 2 , (4.20) 

where 

«1 = EE e' adPoot {a,d) , (4.21) 

ad t=p,n 

with 

e' ad = -\T, T,(-l) J+Ja+lb (2J + l)-^=Vj N T=1 (ab, bd)^l + 5 ab )(l + 5 cd ) , (4.22) 

and 

^2 = ^E E ^T(ac,6rf)[^ r (0x/) XT (j)] J=0 . (4.23) 

A abcdK,T=0,l 

Here, we define Pkmt as 

Aft-AfT = PKMp+ {-^fpKMn , (4.24) 

and the E KT are given by 

£k r= o(ac,M) = (_l)i6+ic Ej (_i)J(2j+i)|^ £ ^} ^(1 + ^(1 + ^) 

x I fe =1 (a&, erf) + |(V/ T=0 (a6, erf) - V& =1 (a&, erf))] , (4.25) 



E KT=1 (ac,bd) = _(-i)i»+ic Ej (_i)J(2j + 1 



x 



jd jc # 

I (Vj^ =0 (a&, erf) - Kf T=1 (afo, erf 



[1 + 5 ab )(l + 5 cd ) 



(4.26) 



In this isospin formalism, since Aj T -MTz(ab) = (— 1)>+^ J+T A JT -MTz{ba), the definitions 
of the symmetric and antisymmetric parts of Vj^(ab, cd), Vf T (ab, cd), and Vj^(ab, cd) become 



VjT A (ab, cd) = - 



Vf T (ab, cd) ± (-l) J+ ^ + ^ +T -Vj^(6a,crf) 



(4.27) 
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Note that these expressions allow less freedom in manipulating the decomposition since 
we have to couple proton with proton and neutron with neutron in forming the density 
operators. Also note that ExT=o{0'C,bd) — Ej(T=i{ ac ) bd) is an invariant related only to the 
physical part of the interactions, (Vj^ =1 + V^ =0 ). We can choose all E KT=1 to be zero in 
the above (by setting Vf T=1 = VJ^ =Q ) leaving E KT=0 completely determined by the physical 
matrix elements. In that case, we can halve the number of fields to be integrated. 

If we now diagonalize the E KT (i,j) as before and form the operators 

Qkmt{o) = 7= ■ — ApKMT(a) + {-l) M Pk-mt{u)) , (4.28a ) 

V 2(l+0M0j 

P KMT (a) = - * (pKM T {a) - (-l) M p K -M T {a)) , (4.28b ) 

the two-body part of the Hamiltonian can finally be written as 

H2 ^EEM") E (Q 2 KMT^) + P 2 KMT {a)) . (4.29) 

Z KT a M>0 

In this decomposition, the one-body hamiltonian h of the HS transformation does not mix 
protons and neutrons. We can then represent the proton and neutron wavefunctions by 
separate determinants, and the number of neutrons and protons will be conserved rigorously 
during each Monte Carlo sample. For general interactions, even if we choose nonzero Ekt=i 
matrix elements, the number of fields involved is half that for a decomposition that mixes 
neutrons and protons, and the matrix dimension is also halved. These two factors combine to 
speed up the computation significantly. In this sense, an isospin formalism is more favorable, 
although at the cost of limiting the degrees of freedom embodied in the symmetric matrix 
elements Vj. 



5 Monte Carlo quadrature 

The manipulations of the previous sections have reduced the shell model to quadrature. 
That is, thermodynamic expectation values are given as the ratio of two multidimensional 
integrals over the auxiliary fields. The dimension D of these integrals is of order N^N t , which 
can exceed 10 5 for the problems of interest. Monte Carlo methods are the only practical 
means of evaluating such integrals. In this section, we review those aspects of Monte Carlo 



quadrature relevant to the task at hand. Our discussion is adapted from that of Ref. |p0 
where a more general exposition can be found. 

Monte Carlo quadrature can be a very efficient way of evaluating integrals of high dimen- 
sion. The name "Monte Carlo" arises from the random or "chance" character of the method 
and the famous casino in Monaco. The essential idea is not to evaluate the integrand at ev- 
ery one of a large number of quadrature points, but rather at only a representative random 
sampling of fields. This is analogous to predicting the results of an election on the basis of 
a poll of a small number of voters. 

To apply Monte Carlo quadrature to the shell model (i.e., expressions such as (3.19)), 
the weight function W a must be real and non-negative. While the properties of W clearly 
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depend upon the interaction and decomposition used, these conditions are not satisfied 
for the natural decompositions of most realistic hamiltonians. However, a broad class of 
schematic interactions and all realistic interactions satisfy these conditions closely enough 
(or, occasionally, exactly) so that Monte Carlo quadrature can be used. Thus, for the 
purposes of the present discussion, we assume W a > and postpone a discussion of the 
so-called "sign problem" to section 7 below. 
To understand the Monte Carlo method, we recast the ratio of integrals in Eq. (3.19) as 

(Cl) = J d D aP a n a , (5.1) 

where 

" W " (5.2) 



Jd D aW a ' 

Since / d D aP a = 1 and P a > 0, we can think of P a as a probability density and (Cl) as the 
average of fl a weighted by P a . Thus, if {a s , s — 1, . . . , S} are a set of S field configurations 
randomly chosen with probability density P a , we can approximate (Cl) as 

(n)«i£ n " (5.3) 

D 8=1 

where fl s is the value of fl a at the field configuration a s . Since this estimate of (Cl) depends 
upon the randomly chosen field configurations, it too will be a random variable whose average 
value is the required integral. To quantify the uncertainty of this estimate, we consider each 
of the fh as a random variable and invoke the central limit theorem to obtain 



2 - L Jd D aP a (Cl (T -(Cl)) 2 ^^J2(n s -(Cl)) 2 . (5.4) 



s 



<«> " S 



Equation (5.4) reveals two very important aspects of Monte Carlo quadrature. First, the 
uncertainty in the estimate of the integral decreases as S^ 1 ^ 2 . Hence, if more samples are 
used, we will get a more precise answer, although the error decreases very slowly with the 
number of samples (a factor of four more numerical work is required to halve the uncertainty). 
The second important point to realize from (5.4) is that the precision of the result depends 
upon the extent to which Q a deviates from its average value over the region of integration; 
the best case is when Q a is as smooth as possible. 

The discussion above shows that Monte Carlo quadrature involves two basic operations: 
generating field configurations randomly distributed according to the specified distribution 
P a (or, equivalently, W a , to which it is proportional) and then evaluating the observable 
Q a for these fields. The second operation is straightforward (details are discussed in the 
following section and the appendix), but it is less obvious how to generate the required field 
configurations. 

When the probability distribution has a simple analytical form, there are a number of 
methods that can be invoked to generate independent random samples [0 . However, in the 



present case P a has a complicated dependence upon the auxiliary fields that is not amenable 
to these methods. Rather, we employ the Metropolis, Rosenbluth, Rosenbluth, Teller, and 
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Teller algorithm [fjEfl , which requires only the ability to calculate the weight function for a 
given value of the integration variables. 

Although the algorithm of Metropolis et al. can be implemented in a variety of ways, we 
begin by describing one simple realization. Suppose that we want to generate a set of samples 
distributed according to a (not necessarily normalized) weight function W a . The Metropolis 
algorithm generates an arbitrarily long sequence of samples {a k ,k = 0,1,2,...} as those 
field configurations visited successively by a random walker moving through cr-space; as the 
walk becomes longer and longer, the points it connects approximate more closely the desired 
distribution. 

The rules by which the random walk proceeds through cr-space are as follows. Suppose 
that the walker is at a point a k in the sequence. To generate cr fc+1 , it makes a trial step 
to a new point Ot- This new point can be chosen in any convenient manner, for example 
uniformly at random within a multi-dimensional cube of small side 5 about o k . This trial 
step is then "accepted" or "rejected" according to the ratio 

r = w t • (5 - 5 > 

where W t is W evaluated at a t and W k is W evaluated at a k . If r is larger than one, then the 
step is "accepted" (i.e., we put a k+ i = at), while if r is less than one, the step is accepted 
with probability r. This latter is conveniently accomplished by comparing r with a random 
number r\ uniformly distributed in the interval [0, 1] and accepting the step if 77 < r. If the 
trial step is not accepted, then it is "rejected," and we put a k+ i = cr k . These rules generate 
0fc + i, and we can then go on to generate a k +2 by the same process, making a trial step from 
Ofc+i- Any arbitrary point, cr , can be used as the starting point for the random walk. 

To prove that the algorithm described above does indeed generate a sequence of points 
distributed according to W, let us consider a large number of walkers starting from different 
initial points and moving independently through cr-space. If N k (a) is the density of these 
walkers at step k, then the net number of walkers moving from field configuration a to field 
configuration p in the next step is 



AN(a) = N k (a)P(a - p) - N k (p)P(p -> a 

~ N k (a) _ P(p^a) 
N k (p) P{o^ P ) 



N k (p)P(a -+p)\m- ^} • (5.6) 



Here, P(a —* p) is the probability that a walker will make a transition to p if it is at a. This 
equation shows that there is equilibrium (no net change in population) when 

N k (<r) P(p^cr) JUa) 

N k (p) P(a - p) N e (p) 1 • ) 

and that changes in N{a) when the system is not in equilibrium tend to drive it toward 
equilibrium (i.e., AA^(cr) is positive if there are "too many" walkers at a, or if Nk{a)/N k {p) 
is greater than its equilibrium value). Hence it is plausible (and can be proved) that, after 
a large number of steps, the population of the walkers will settle down to its equilibrium 
distribution, N e . 
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It remains to be shown that the transition probabilities of the Metropolis algorithm lead 
to an equilibrium distribution of walkers N e (a) ~ W a . For the random walk governed by 
the rules described above, the probability of making a step from cr to p is 

P(a - p) = T(a - p)A(a - p) , (5.8) 

where T is the probability of making a trial step from cr to p and A is the probability of 
accepting that step. If p can be reached from a in a single step (i.e., if it is within a cube of 
side 5 centered about cr), then T(cr — > p) — T(p — > cr), so that the equilibrium distribution 
of the Metropolis random walkers satisfies 

N e (a) _ A(p -> cr) 

7V e (p) A(a - p) • l °' yJ 

If W ff > W^, then A(p -> a) = 1 and A(a -> p) = H/p/Wx, while if W ff < W 7 , then 
A(p — > er) = W a /W p and A(cr — > p) = 1. Hence, in either case, the equilibrium population 
of Metropolis walkers satisfies 

W - ^ (5 10) 

N e (p) ~ W p (5 - 10) 

so that the walkers are indeed distributed with the correct distribution. 

Note that although we made the discussion concrete by choosing cr t in the neighborhood 
of cr fc , we can use any transition and acceptance rules that satisfy 

T(p -> a)A(p -> _W a 

T(a -> p)A(cr -> p) W p " 1 ' ' 

Indeed, one limiting choice is T(a — > p) = W^, independent of cr, and A = 1. This is the 
most efficient choice, as no trial steps are "wasted" through rejection. However, this choice 
is somewhat impractical, because if we knew how to sample W to take the trial step, we 
wouldn't need to use the algorithm to begin with. 

An obvious question is "If trial steps are to be taken within a neighborhood of cr*., how do 
we choose the step size, 5?" To answer this, suppose that cr^ is at a maximum of W, the 
most likely place for it to be. If 5 is large, then W t will likely be very much smaller than 
W k and most trial steps will be rejected, leading to an inefficient sampling of W . If 5 is very 
small, most trial steps will be accepted, but the random walker will never move very far, and 
so also lead to a poor sampling of the distribution. A good rule of thumb is that the size of 
the trial step should be chosen so that about half of the trial steps are accepted. 

One bane of applying the Metropolis algorithm is that the field configurations that make 
up the random walk, a ,a±, . . ., are not independent of one another, simply from the way 
in which they were generated; that is, (Jk+i is likely to be in the neighborhood of a k . Thus, 
while the samples might be distributed properly as the walk becomes very long, they are 
not statistically independent of one another, and some care must be taken in using them to 
evaluate observables. In particular, if we estimate (Cl) from Eq. (5.3) using the successive 
field configurations in the random walk, the estimate of the variance given by Eq. (5.4) is 
invalid because the Q k are not statistically independent. This error can be quantified by 
calculating the auto-correlation function 

_ (n k n k+i ) - w 2 r _ 19 . 
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where (. . .) indicates average over the random walk. Of course, Co = 1, but the non-vanishing 
of Cj for % > means that the fVs are not independent. What can be done in practice is to 
compute the observable and its variance using samples along the random walk separated by 
a fixed interval, the interval being chosen so that there is effectively no correlation between 
the field configurations used. An appropriate sampling interval can be estimated from the 
value of % for which C becomes small (say < 0.1). 

Another issue in applying the Metropolis algorithm is where to start the random walk; i.e., 
what to take for uq. In principle, any location is suitable and the results will be independent 
of this choice, as the walker will "thermalize" after some number of steps. In practice, an 
appropriate starting point is a probable one, where W is large. Some number of thermal- 
ization steps then can be taken before actual sampling begins to remove any dependence on 
the starting point. 

6 Numerical methods and computational issues 
6.1 Number projection 

For large quantum systems in which many valence fermions participate, the grand canonical 
ensemble (with a fixed chemical potential, but fluctuating particle number) is perfectly ac- 
ceptable, and simple formulas exist for calculating the required one-body observables (e.g., 
Eqs. (A. 15,16)). But in nuclei, where only a few valence fermions are important in deter- 
mining the low lying states, and the properties of neighboring even- A and odd- A systems 
differ dramatically, the canonical (fixed-number) ensemble must be used. In this section, 
we discuss how to project the fixed-number trace from the grand-canonical ensemble. For 
simplicity, we will consider only a single species of fermions with number operator N; the 
required generalization to two species (protons and neutrons) is straightforward, as the HS 
transformation used does not involve density operators that mix the species, so that U a 
factors into separate evolution operators for neutrons and protons. 

The grand canonical partition function for a one-body evolution operator U a is given by 
Eq. (A.12) as 

Z = Tre^U a = det(l + e^U a ) , (6.1) 

where \i is the chemical potential. Note that since the matrix corresponding to N is simply 
the unit matrix, we may include the chemical potential term in each of the h defining U a , or 
write it as a common prefactor, as we have done in (6.1). If we write the N s eigenvalues of 
Uo- as e~@ £x , where e\ may be complex, we can express (6.1) as 

Z = H(l + e^- £ ^) , (6.2) 

A 

which is a form reminiscent of the usual fermi-gas expressions. 
The corresponding canonical partition function for A particles in the valence space is 

Z A = Tr A U a = TrP A U a , (6.3) 

where P A = 5(A — N) is the number projector. To effect the operation of P A) we can write it 
as the exponential of the one-body operator N using the fact that N has integer eigenvalues 
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0,1,. ..,A S : 

P A = e -^ A ^ e -iM e (^+i4>)N ^ (6 4) 

Jo 2tx 

where ji is an arbitrary c-number, which will be specified below. 
The canonical trace thus becomes 

Z A = e~^ A r^e-^Tre^+^Ur (6.5a) 
Jo 2ix 

= e -/3M [ 27T TT (1 + e i^-e x)) (6>5b ) 

Jo 2ir > 



From (6.5b), we can see that the effect of the number projection is to sum all different 
products of A terms of the form e~ l3ex . Similar expressions can be written for expectation 
values of few-body operators following Eqs. (A. 15, A. 16). For example, for a one-body 
operator Cl, 

Tr A U a Q = e~^ A ^ det(l + e^UJtr } e^U^ . (6.6) 

The parameter \i is arbitrary and hence can be chosen for numerical convenience. Since 
the grand-canonical trace in the integrand of (6.5b) contains contributions from canonical 
ensembles with all < A < N s , and these vary strongly with A, fi should be chosen 
to emphasize that contribution from the particular value of A desired; a poor choice can 
lead to numerical difficulties in the projection. This can be done in analogy with usual 
thermodynamic treatment. In particular, if the eigenvalues of U CT are ordered so that Re^i < 
Ree 2 < . . . < Ree Ns , then a good choice for /i is (Ree^ + Re£A+i)/2. 

The numerical effort implied by Eqs. (6.5) or (6.6) need not be intimidating. First, because 
the integer eigenvalues of A range from to N s , the maximal Fourier component in the <fi 
integration is e lNs ^, so that an Appoint quadrature can be exact. That is, for /(0) one of 
the required integrands, 



L 



* m=l * 



Second, the matrix algebra associated with expressions like (6.6) is minimized by computing 
only once both the eigenvalues e~ l3ex and the diagonalizing transformation matrix T a \ of 
Uo- and then using these to evaluate the integrand at each quadrature point <f) m . Thus, for 
example, (6.5b) is simply evaluated, and the matrix trace required in the integral of (6.6) is 
evaluated as 



£ l + e ^-e^ ^ £x)+l " £ t^cwt^ . (6. 

A=l 1 ^ C a/3=l 



6.2 Auxiliary field quadrature 

When the number of auxiliary field variables becomes large, naive application of the sampling 
algorithm of Metropolis et al. as described in Section 5 becomes inefficient. For example, in 
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the 170 Dy calculation described in Section 8.8 where there are some 10 5 fields, when the step 
size S (taken to be the same for each field variable) is adjusted so that the acceptance ratio is 
approximately 0.5, the auto correlation length for the energy (as computed from Eq. (5.12)) 
is more than 200 sweeps (i.e., more than 200 updates of the field variables at each time slice). 

The Metropolis efficiency in generating uncorrelated field configurations can be improved 
significantly by approximating the continuous integral over each a an by a discrete sum derived 



from a Gaussian quadrature |59]. In particular, the relation 

J —oo 

is satisfied through terms in (A/3) 2 if 

f(a) = \ [8{a - a ) + 8{a + <r ) + 45(a)] , (6.10) 
6 

where o"o = (3/VA/3) 1 / 2 . (Note that the commutator terms render the HS transformation 
accurate only through order A/3 anyway.) 

In this way, each a an becomes a 3-state variable and the integrals in (3.19) become (very 
large) sums, which can be sampled using the algorithm of Metropolis et al. In particular, 
to update the fields at a given time slice, we select at random a fraction q of them and 
assign each of these to (— cr , 0, or cr ) with probability (1/6, 2/3, and 1/6). To satisfy the 
condition (5.11), the move is then accepted or rejected according to the ratio of the old and 
trial values of TrxU a . The fraction q is adjusted to give an acceptance ratio of about 0.5. 
When used in the 170 Dy calculation, this discretization reduces the energy autocorrelation 
length to only five sweeps, a factor of 40 improvement relative to the continuous case, with 
no loss of accuracy. Similar, albeit less substantial, savings accrue in smaller calculations. 

6.3 Stabilization at low temperatures 

The excitation energy of the first excited state in even-even nuclei with A < 80 is 1-2 MeV. 
Thus, one may calculate at a temperature of T = 0.5 MeV to get a reasonable description of 
the ground state. However, in odd- A and odd-odd systems the excitation energy can be as 
low as tens of keV, so that cooling to the ground state requires calculations at /3 > 2.0 MeV -1 . 
Unfortunately, numerical instabilities arise for larger (3 due to the multiplication of many 
U matrices together; the product becomes increasingly illconditioned, with exponentially 
divergent numerical scales, as illustrated below. 

Physically, the one-body evolution operator U n amplifies low-energy states for any par- 
ticular field configuration, and attenuates high-energy states. The most important states, 
those near some intermediate Fermi surface, are buried exponentially by the states at the 
bottom of the spectrum, and extraction of relavent information becomes increasingly dif- 
ficult with lower temperatures. This technical problem, resulting from the limits of finite 
numerical precision achievable on any digital computer, can be circumvented using singular 
value decompositon (SVD) methods, which were first applied in Green's Function Monte 
Carlo simulations of interacting electron systems 5"3fl . Here we adapt the SVD methods to 
the canonical calculations of the nuclear problem. 
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Let us suppose that a matrix A can be decomposed into three matrices, A = SDV, where 
S^S = 1, V^V = 1, and D is a real, diagonal matrix. Schematically, this decomposition 
looks like 
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(6.11) 



where the size of the symbols (X, x, x) indicates the magnitude of the elements. Note that S 
and V are well controlled matrices shown schematically as having unit scale, and that all the 
scales of the original matrix A are contained in D. Note also that once the multiplication 
has been performed all elements of the resulting matrix are of the same (largest) scale, and 
the product is essentially an outer product of the first column of S and the first row of V. 
The smallest numerical scales exist only implicitly as differences of large matrix elements, 
and as such are difficult to recover on a computer with finite precision. 

In contrast, a matrix that explicitly displays its small scales, such as a column stratified 
matrix M, can be factored stably. 
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Multiplication on the left of M by a transformation matrix combines only elements in a 
given column which are all of the same scale. Thus there is no loss of information. Multi- 
plication on the right by V -1 combines columns of different scales, but does not overwrite 
any small scale as long as large scaled columns are scaled down appropriately before the 
addition. 

These problems are quantified by the condition number, defined as the ratio of the largest 
element of D to the smallest after the transformation has taken place. When this ratio 
is greater than the machine accuracy (usually about 10 12 for double precision), then the 
calculation will become unstable, and techniques of stable matrix multiplication must be 
employed. Numerical problems first occur in the multiplication of U„ matrices to create 
U = n^ r i 1 U„. Although each U n has a good condition number, the final product can have 
an unacceptable condition number. The above examples can be combined to disentangle the 
scales inherent in the multiplication of the U n matrices. 

SVD matrix multiplication works by isolating the scales of the problem at each multipli- 
cation step. Let the matrix W n = n^ =1 *U m , be the partial product of the matrix multipli- 
cations, with W = 1, Wi = Ui, W 2 = U 2 Ui, etc. Suppose we have already decomposed 
the matrix W n in SVD format. We want to multiply the well-conditioned matrix U n with 
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W n _i. The multiplication proceeds as follows: 











U n W n _i = (U n S n _iD n _i) V n _i = 


U n S n _i 















Vn-l ■ (6.13) 



Multiplication of the matrices in brackets yields a column stratified matrix that can be 
put into SVD format (i.e., S'D'V) without loss of scale. Thus 



(X 



U n W n _! = 



X 



X x x 
\X x x J 



V n _, = (S'D'VOV^! = S'D'tV'V^) = S n D n V n , (6.14) 



where in the last step we have defined S n = S', D n = D', and V n = V'V n _i. Provided 
that the V matrices are sufficiently well-conditioned that we can multiply many of them 
together without scale problems, we are able to perform stable matrix multiplication of all 
the one-body evolution operators. 

For evaluation of the quantum mechanical trace, we need to calculate 1 + XJ a consistent 
with SVD. We assume that SVD has previously been performed on U CT such that = SDV. 
The goal is to keep all of the large scales in the problem separated from the small ones. This 
is done by calculating 

det (1 + U CT ) = det [S (S^V- 1 + D) v] = det SSDVV = det S' det D' det V . (6.15) 

Here we have performed a new SVD on the matrix in parenthesis to get SDV (after the 
addition of S~ 1 V~ 1 to D), and we define S' = SS and V = VV The other quantity we 
need is 



;i + u tr )- 1 u ff = i- 



1 _ v'^D'^S 



i-i 



(6.16) 



1 + U CT 

where we have used the inverses of the matrices calculated in Eq. (6.15). These steps can be 
carried over exactly into the number projection method described in section (6.1). 

As an example we consider the matrix multiplications to obtain U CT for 54 Fe in the course 
of the calculations described in section 8.1. Shown in Fig. 6.1 is an example of the condition 
number with and without SVD for various (3, with Af3 = 1/16 MeV -1 . Numerical round-off 
causes the condition number to virtually saturate in the non-SVD case. Shown in Fig. 6.2 
are the elements of the diagonal matrix D using SVD for all multiplies, compared to using 
SVD only on the full U CT . Note that the symmetry of U CT (that eigenvalues come in complex 
conjugate pairs, see section 7.1) is evident when the matrix-multiplies are accomplished 
through SVD, where the real elements of D come in pairs. This is not the case without SVD 
multiplication. The numerical eigenvalues of violate the symmetry only after (3 > 2.0 
(and at larger f3 in larger systems). 

We have made limited use of SVD in the calculations presented in this Report. In test 
calculations we find that we can increase (3 almost indefinitely (e.g. to (3 = 6 in 22 Ne) with 
SVD, whereas without the stabilization, the maximum attainable (3 is 2.0 MeV -1 . In the 
170 Dy and Xe calculations described in section 8, f3 — 5 MeV -1 is attainable without SVD. 
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6.4 Maximum entropy techniques 



The response functions that can be evaluated in SMMC are given by Eq. (3.6). Calculation 
of the strength function Sq(E) for a given operator can be difficult due to the ill-posed 
inversion of the Laplace transform required by Eq. (3.8a). In this section, we discuss one 
particular inversion technique, Maximum Entropy (MaxEnt), which we have used to find 
the Gamow- Teller strength functions presented in section 8.6. 
We wish to find the function S(E) from the function R(t) using the kernel K(t, E): 



R(t) = J dEK{j, E)S(E) . 



(6.17) 



The function R(t) is known at N grid points Tj (i = 1, ■ ■ • , N), where the values of R are 
given at T{ by r.; with errors cr.j. Fluctuations in different data points may be correlated, 
in which case the covariance matrix Cy = ((r^ — fj)(rj — f,)) is diagonalized to find the 
independent modes, and and the kernel are then transformed to this new basis. The goal 
is to find the best S(E) by analysing the \ 2 figure of merit, 



x 2 {s} = E 



'n-Riisy 



<7a 



(6.18) 



of the data r; from the fit values Ri produced by the trial inverse in the points Tj. 



Ri{S} = j dEK{r h E)S{E). 



(6.19) 



Direct minimization of x 2 is numerically stable in only the simplest of circumstances. Most 
of the modern ways to regularize the problem are based on combining \ 2 with some other 
auxiliary well-conditioned functional, P{S}. This functional should have a minimum at 
a smooth function S(E) and assign a high cost to strongly oscillating functions. Thus we 
minimize the joint functional 

\x 2 {S} + P{S} . (6.20) 
We use the information-theoretic entropy to define P {S} such that 



P{S} = a J 



dE 



m(E) - S(E) + S(E) In 



( S(E) \ 
\m(E)J 



(6.21) 



where m(E), the default model, and a are chosen to describe the a priori knowledge about the 
shape of the original S(E). The choice of m(E) and a is more doubtful than the MaxEnt 
technique itself. In our application to the Gamow- Teller strength functions we make the 
simplifying assumption that m(E) is a Gaussian function with appropriately chosen center 
(see section 8.6), and a — [J dEm(E)}^ 1 . 

In order to minimize the functional (6.20) we employ the technique of Ref. p4]| , which 
involves an iterative sequence of linear programming problems. To do this we first expand 
Eq. (6.21) to second order in S(E) about some positive function f(E), to obtain 



P{f\S} = aJ 



dE 



m — 



+ 




- 1 




(6.22) 
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If the true minimum S(E) of the non-quadratic functional in Eq. (6.21) is taken as a point 
of expansion f(E) in (6.22), then it also gives the minimum of the corresponding quadratic 
functional 



S(E) 



mm 



- x *{a} + P{S\a} 



(6.23) 



We transform this identity to a recursive equation that, upon iteration, leads to the desired 
solution. To ensure positivity we keep a fraction of the result of the previous iteration in the 
next one, and minimize the functional with the restriction that S > 0. Defining n as the 
iteration step, we obtain, 



mm 

5>0 



1 



{S} + p{fV> I s} 



with 



f^(E) = ^S ( - n - 1 \E) + (l-^ n \E), 
and the default model is taken as the starting approximation to S, 

S ( - 0) (E) = S ( - 1) (E)=m(E) . 



(6.24) 
(6.25) 
(6.26) 



The rate of convergence and stability are controlled by the mixing parameter £; a value of 
£ = 0.3 is a reasonable choice to gaurantee stability. 

Each iteration is a linear regularization task with a quadratic regularization functional 
reinforced by an additional stabilizing effect of the constraints S(E) > 0. Thus, each iteration 
gives a nonnegative function S^ n '(E) and a positive definite function f( n \E). Test cases of 
simulated data for which the answer is known give convergence within twenty iterations. 

We have used the technique described above to obtain the strength functions from the 
Gamow- Teller response functions, as presented in section 8.6. A different approach to Max- 



Ent was outlined in p7| , but has proven less useful due to difficulties in obtaining a converged 
value of a. The advantage of the present method is that knowledge of a is not necessary, 
and that the scheme provides an iterative method of converging to the most likely solution. 



6.5 Computational considerations 

The present SMMC code is roughly the fourth major revision of a program whose develop- 
ment began in late 1990 with a single-species, single-j-shell version. It is now a modular 
package of some 10,000 commented lines of FORTRAN. All floating point computations 
are double precision (64-bit). 

The package performs all of the functions necessary for shell model Monte Carlo cal- 
culations: initialization, thermalization of the Metropolis walk, generation of the Monte 
Carlo samples, evaluation of static observables and response functions (canonical or grand- 
canonical), MaxEnt extraction of strength functions, and the extrapolation in g required 
to solve the sign problem, as discussed in section 7. Samples may be analyzed "on-line" 
and/or stored for post-processing in subsequent runs. The data input and results output use 
standard shell model conventions, so that it is easy to change the two-body matrix elements 
of the interaction, to incorporate additional one- or two-body observables in the analysis, or 
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to add or change the orbitals in the calculation. The code has been de-bugged and tested 
extensively against direct diagonalization results in the sd- and lower ^/-shells. Its opera- 
tion by an experienced user can be described as "routine," although it takes several weeks 
to acquire that experience. 

The SMMC package is highly portable. Initial development was largely under VAX/VMS. 
Subsequently, we have ported the code to DEC ALPHA'S and HP 730's operating under 
Unix, the Intel i860 CPU used in the DELTA and PARAGON parallel machines, the IBM 
RS-6000 processor used in both workstations and the SP-1 and SP-2 parallel machines, and 
the Fujitsu VPP500 shared memory machine (at RIKEN). Very few problems are encountered 
in porting to a new machine, and the operation generally takes less than a day. 

Shell model Monte Carlo calculations are extraordinarily well-suited to Multiple Instruc- 
tion/Multiple Data (MIMD) architectures. Indeed, our code is "embarrassingly parallel": 
separate Metropolis random walks are started on each computational node, which then pro- 
duces a specified number of Monte Carlo samples at regular intervals during the walk. Data 
from all of the nodes are sent to a central node for evaluation of the Monte Carlo averages 
and their uncertainties, and perhaps for storage of the sampled field configurations in a file 
for post- analysis. The post-processing of the stored samples is also divided among the nodes. 

To date, we have implemented the parallel version of our code on the Intel DELTA, and 
PARAGON machines at Caltech (each with 512 i860 processors), on the 128-processor IBM 
SP-1 at ANL, the 512-processor IBM SP-2 at Maui, and on a Fujitsu VPP500 shared memory 
vector processor (24 CPU's). In all cases, the ratio of communications to computation is 
very low, with efficiencies always greater than 95%. 

To circumvent the limited memory on the DELTA (12.5 MBytes) nodes, we have also 
produced a version of the code in which the chains of U a matrices effecting the time evolution 
are split over two or more nodes, information being passed between them only at the time 
slice where they "join." However, such an implementation is a significant sacrifice in speed, 
as while one processor updates a chain, the other processors working on the same chain must 
sit idle. 

Table 1 shows benchmarks of our code on various single processors. The test calculation 
involved a canonical ensemble in the full p /-shell ( N s = 20 for each type of nucleon, implying 
20 x 20 matrices) using a realistic interaction. N t = 32 time slices were used at (5 — 2 MeV -1 
(A/3 = 0.0625 MeV -1 ). Thirty static observables and seven dynamical response functions 
were calculated at a single g- value (see section 7). Note that the computational speed is 
independent of the interaction and of the number of nucleons occupying the shell. Beyond 
using library subroutines (BLAS and LAPACK routines), no attempt was made to optimize 
the assembly level code in any of these cases. Approximately 40% of the computational effort 
is in matrix-multiplies. A significant remainder of the effort goes into building the one-body 
Hamiltonian (13%), setting up the one-body evolution operators (15%), and calculating 
two-body observables (15%), none of which is easily vectorizable. 

In general, the computation time scales as N^N t , and is spent roughly equally on the 
dynamical response functions and the static observable sampling. Of course, the stabilization 
discussed above increases these times significantly (by at least a factor of 5). 

The memory required for our calculations scales as N^N t . is the average of the squares 
of the numbers of neutron and proton single-particle states. Sample values of N s for one 
isospin type are shown for various model spaces in Table 2. 
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Processor 


Peak MF 


Average MF 


Samples/hr. 


i860 


35 


9 


44 


IBM-SP2 Thin66 


56 


36 


179 


ALPHA-400 


56 


28 


141 



Table 1: Benchmarks of the SMMC code in various processors. 



Model Space 


N s 


Op 


6 


ls-Od 


12 


lp-Of 


20 


lp-0/-#9/2 


30 


2s-ld-0g 


30 


2p-lf-0g 9/2 -0i 13 /2 


44 



Table 2: Matrix dimension for various model spaces. 



The code is currently structured so that a calculation with N s = 32, six j-orbitals, and 
N t = 64 time slices will fit in 12 MB of memory (note that all field variables are stored as 
1-byte integers because of the 3-point Gauss-Hermite quadrature described in section 6.2). 
A calculation in the (lp-0f)-(2s-ld-0g) basis has N s = 50 and would require about 64 MB 
of memory for Nt = 64 time slices and about 128 MB for Nt = 128. For calculations with a 
"good" Hamiltonian, as defined in the next section, these storage estimates can be reduced 
by a factor of two by exploiting the time reversal symmetry of the U a , albeit at some loss in 
speed. 



6.6 Contrast with interacting electron systems 

It is interesting to contrast the SMMC simulations of nuclear systems with similar work 
on the Hubbard model in the context of high-temperature superconductivity. The latter 
involve electrons hopping locally on a two-dimensional spatial lattice interacting with each 



other through an on-site repulsion. For a review of this field see |24 . 

One significant difference is that physically interesting SMMC calculations can be per- 
formed in smaller model spaces. The N s values quoted in the previous section are signifi- 
cantly smaller than the number of sites needed to approximate the translational invariance of 
the CuO planes. However, the matrices representing the fluctuating one-body hamiltonians 
h in SMMC calculations are dense, as the nuclear hamiltonian does not have the local spatial 
structure of the Hubbard model that gives rise to sparse h matrices. 

Another difference is that relatively higher temperatures are of interest in the nuclear case, 
allowing smaller Nt. For a nominal temperature of T = 300 keV and strength of the two- 
body interaction V = 3 MeV, we have T/V = 0.1. In the high-T c case, T « 10 2 K « 10~ 2 
eV, while V ~ 1 eV, so T/V ~ 10~ 2 , an order of magnitude smaller. 
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FIG 6.1 The condition number of U CT with and without SVD for various (5 is plotted for 
54 Fe. 

The SMMC calculations must be performed in the canonical ensemble, as the effects of 
finite particle number are an important aspect of nuclear structure. This necessitates the 
rather cumbersome number projection described in section 6.1. In contrast, the Hubbard 
simulations can be done in the grand-canonical ensemble. 

A final and most crucial difference between the SMMC and Hubbard hamiltonians is 
that the nuclear interaction is predominantly attractive. This significantly diminishes (and 
in some cases eliminates) the sign problems that plague the condensed matter work. These 
problems, and their resolution in SMMC calculations, are the subject of the following section. 

7 Sign problems 

Although we have briefly alluded to "sign problems" in some of the previous chapters, vir- 
tually all of our discussion has been based on the premise that the weight function W a of 
Eq. (3.19) is non-negative for all field configurations a. Were this not the case, the function 
P a of Eq. (5.2) could not be interpreted as a probability density, and the MC quadrature 
described in section 5 could not be effected. 

Unfortunately, many of the hamiltonians of physical interest suffer from a sign problem, 
in that W a is negative over significant fractions of the integration volume. To understand 
the implications of this, let us rewrite Eq. (5.1) as 

(Cl) = J d D aP a $M* , (7.1) 

where 





w a | 




fd D a 


1 w a 


1 $ 

1 ^cr 
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and $ CT = W a / \ W a | is the sign of the real part of W a . (Note that since the partition function 
is real, we can neglect the imaginary part.) Since | W a | is non-negative by definition, we 
can interpret it, suitably normalized, as a probability density, so that upon rewriting (7.1) 

as 

,K_ fda\w tr \$ n a _ ($n) r79 , 
fd<j\w a \<s> a - (*> ' {7 - Z) 

we can think of the observable as a ratio in which the numerator and denominator can be 
separately evaluated by MC quadrature. Leaving aside the issue of correlations between 
estimates of these two quantities (they can always be evaluated using separate Metropolis 
walks), the fractional variance of (Cl) will be 



<""> ■ 1 2, (7.3) 



which becomes unacceptably large as the average sign (<&) approaches zero. The average 
sign of the weight thus determines the feasibility of naive MC quadrature. Only a handful 
of interacting electron systems are known to give rise to a positive-definite path integral 



2% : the one-dimensional Hubbard model, the half-filled Hubbard model, and the attractive 
Hubbard model at any dimension and filling. 

When (<&) is near unity, we can attempt MC evaluation of the numerator and denominator 
separately, as will be illustrated in the example of section 8.8 below. However, for the most 
general case, we must have a deeper understanding of the sign problem. We rewrite the 
canonical expectation value of an observable fl as 

c Tr(Cle-P 6 ) ^ [DaWMh 

{ ' Tr(e-f> 6 ) ~ JDaW a <S> a ' 1 ' ' 
where in the spirit of Eq. (7.1,7.2) we have introduced a positive definite weight 

W a = G a \ Tr A f> CT | , (7.5) 



and the Monte Carlo sign 



= J^Sf- . (7.6) 



The sign problem arises because the one-body partition function Tr aUo- is not necessarily 
positive, so that the Monte Carlo uncertainty in the denominator of Eq. (7.4) (the W- 
weighted average sign, ($)) can become comparable to or larger than (<&) itself. In most 



cases (<&) decreases exponentially with (3 or with the number of time slices [53 



7.1 Hamiltonians without sign problems 

An important class of interactions (pairing+quadrupole) free from the sign problem (i.e., 
$0- = 1) was found in Ref. ^7|. These are realized when V a < for all a in Eq. (3.16). In 

that case, s a — 1 for all a, and, in the linearized hamiltonian given by Eq. (3.19c), O and O 
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couple to complex conjugate fields. In order to understand the implications of this situation, 
we represent the single-particle wave function by an ^-component vector of the form 



jm 
jfh 



(7.7) 



with Ng/2 states m > in the first half of the vector, and their time reversed orbitals in 
the second half. Due to Eq. (3.16) and the fact that time-reversed operators are coupled to 
complex-conjugate fields, the matrix h n has the structure 



A R 



-b; a; 

and one can easily verify that the total evolution matrix 



(7.8) 



U CT = IJexp (-KA(3) = ( _q* £ ) , (7.9) 

is of the same form. Here A, B, P, and Q are matrices of dimension N s /2. One can show that 
this matrix has pairs of complex-conjugate eigenvalues (e,e*), with respective eigenvectors 

an d ( ^ )■ When e is real, it is two-fold degenerate, since the two eigenvectors are 

distinct. 

For the grand-canonical ensemble, the grand-canonical trace is given by 



TrU = det 



1 + 



P Q- 

-Q* P* 



U(l+e x )(l + e* x )>0. (7.10) 



A 



If only particle- type conserving (neutron-proton) density operators are present in Eq. (3.16), 
each type of nucleon is represented by a separate Slater determinant having the structure 
(7.9) and therefore TrU p ■ TrU n > since TrU p > and TrU n > 0. 

In the zero-temperature formalism, if the trial wave function for an even number of particles 
is chosen to consist of time-reversed pairs of single particle states, 

where a and b are matrices with dimension x ^j, then ^TJtfr is a N v x N v matrix 
with the structure (7.9), and the trace satisfies 



TrU = det 



$t u$ 



(7.12) 



If only particle-type conserving operators are present, then time-reversed pairs of trial wave 
functions can be chosen for both protons and neutrons in an even-even nucleus, giving rise 
to TrU = TrU p ■ TrU n > 0. For N = Z odd-odd nuclei we can choose the final neutron to be 
in the time reversed orbit of the final proton, and the trace is once again positive. 

Eq. (6.5b) embodies a proof that for even-even nuclei (and N = Z odd-odd systems) the 
partition function is real and positive definite under the condition that all V a < 0. This is 
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most easily seen for the case where N s = 2, and A = 1 or A = 2. From Eq. (6.5b), and using 
the iVs-point quadrature of Eq. (6.7), 



Z A = - ]T e (~^) f[[l + e ^) e (^e. 

2 m =l A=l 
1 2 

= _ | e ( ~ i<?imA) + el-^^-^el-^ 1 ' + e [ili>m( - A ~ 1)] Re e (_/3ei) } , (7.13) 

^ m=l 

where we have ignored the e^ 13 ^ term for simplicity, and <p m = mix. It is easy to show 
that for A = 2 the first and second terms survive the sum and that these terms are positive 
definite. It is also easy to show that for A = 1 only the third term, Re e ( -~ /3£l - ) remains 
after the sum is taken. While this term is real, it is not positive definite. Hence there 
appear two distinct sign problems. For even-even nuclei and all V a < 0, the sign problem 
can be overcome; however, for odd- A nuclei, even with all V a < 0, the sign problem remains 
since the contributions to ($) are either positive (+1) or negative (—1), but not complex. 
This demonstration can be extended to all A in a given space. As in the zero-temperature 
formalism, N = Z odd-odd nuclei pose no sign problem, since the evolution is the same for 
both the protons and the neutrons, and the product of the traces is thus positive definite. 

7.2 Practical solution to the sign problem 

For an arbitrary hamiltonian, we are not guaranteed that all V a < 0. However, we may expect 
that a realistic hamiltonian will be dominated by terms like those of the schematic force (this 
is, after all, why the schematic forces were developed) so that it is, in some sense close to a 
hamiltonian for which the MC is directly applicable. Thus, the "practical solution" to the 



sign problem presented in Ref. [29] is based on an extrapolation of observables calculated for 
a "nearby" family of Hamiltonians whose integrands have a positive sign. Success depends 
crucially upon the degree of extrapolation required. Empirically, one finds that, for all of 
the many realistic interactions tested in the sd- and ^/-shells, the extrapolation required is 
modest, amounting to a factor-of-two variation in the isovector monopole pairing strength. 

Based on the above observation, it is possible to decompose H in (3.16) into its "good" 
and "bad" parts, H = Hq + Hb, with 

h g = ^2(^ + e a o a ) + ^ v a {6 a ,d a } 

a Z V a <0 

H B = ~ £ V a {6 a ,d a }. (7.14) 
1 v a >o 

The "good" Hamiltonian Hq includes, in addition to the one-body terms, all the two-body 
interactions with V a < 0, while the "bad" Hamiltonian Hb contains all interactions with 
V a > 0. By construction, calculations with Hq alone have $ CT = 1 and are thus free of the 
sign problem. 

We define a family of Hamiltonians H g that depend on a continuous real parameter g 
as H g = f(g)Hc + gHs, so that H g= \ = H, and f(g) is a function with /(l) = 1 and 
f(g < 0) > that can be chosen to make the extrapolations less severe. (In practical 
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applications f(g) = 1 — (1 — g)/x with x ~ 4 has been found to be a good choice.) If 
the V a that are large in magnitude are "good," we expect that H g =o = Hq is a reasonable 
starting point for the calculation of an observable (fl). One might then hope to calculate 
(Q) g = Tr (Cle~P 9 )/Tr (e - ^ 9 ) for small g > and then to extrapolate to g = 1, but typically 
($) collapses even for small positive g. However, it is evident from our construction that 
H g is characterized by $ CT = 1 for any g < 0, since all the "bad" V a (> 0) are replaced by 
"good" gV a < 0. We can therefore calculate (Cl) g for any g < by a Monte Carlo sampling 
that is free of the sign problem. If (Cl) g is a smooth function of g, it should then be possible 
to extrapolate to g = 1 (i.e., to the original Hamiltonian) from g < 0. We emphasize that 
g = is not expected to be a singular point of it is special only in the Monte Carlo 

evaluation. 

As described in Section 4, the matrices E\ T are constructed from the two-body matrix 
elements VJ T (ab,cd) of good angular momentum J, isospin T, and parity n = (— l)' a+ ^ 
through a Pandya transformation. For interactions that are time-reversal invariant and 
conserve parity, the E^ T (i,j) (here n = (— l) Io+Jc ) are real symmetric matrices that can be 
diagonalized by a real orthogonal transformation. The eigenvectors Pkm(&) play the role of 
O a in Eq. (3.16), and the eigenvalues A^ 7r (a) are proportional to V a . In the Condon-Shortley 
|Kj convention p KM = 7r(— ) k+m Pk-m so that the "good" eigenvalues satisfy sign [A/^a)] = 
ir(-) K+1 |HJ. 

As an example, we consider the mid-pf shell nucleus 54 Fe using the realistic Kuo-Brown 
KB3 interaction Figure 7.1 (upper) shows the eigenvalues Vk-ku = ^(—) K ^Kn(ci) of the 
KB3 interaction; only about half of the eigenvalues are negative. However, those with the 
largest magnitude are all "good." It is possible to use an inverse Pandya transformation to 
calculate the two-body matrix elements VJ T (ab, cd) for the "good" and "bad" interactions, 
allowing the matrix elements of Hq to be compared in Fig. 7.1 (lower) with those of the full 
interaction. The greatest deviation is for J = 0, T = 1 (the monopole pairing interaction), 
where Hq is about twice as attractive as the physical H. In all other channels, Hq and H 
are quite similar. 

In Fig. 7.2 we exemplify the ^-extrapolation procedure for several observables, calculated 
again for 54 Fe. In all cases we use polynomial extrapolations from negative g-values to the 
physical case, g — 1. The degree of the polynomial is usually chosen to be the smallest 
that yields a \ 2 per degree of freedom less than 1. However, in several studies, like the 
one of the p/-shell nuclei reported in section 8.1, we have conservatively chosen second- 
order polynomials for all extrapolations, although in many cases a first-order polynomial 
already resulted in x 2 - values less than 1. At T = the variational principle requires that 
the expectation value of H has a minimum at g = 1. We have incorporated this fact in 
our extrapolations of ground state energies by using a second-order polynomial with zero- 
derivative at g = 1. 

7.3 Validation 

Confidence in SMMC results requires numerical validation of the methods used. Separate 
issues are the validation of the SMMC formalism and algorithms, as developed in sections 
3-6, and of the ^-extrapolation required in calculations with realistic residual interactions 
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FIG 7.1. Decomposition of the Kuo-Brown KB3 interaction. The upper panel shows the 
eigenvalues V a for the various particle-hole angular momenta K. The lower panels show the 
particle-particle two-body matrix elements in selected JT channels. Solid symbols are the 
full hamiltonian, while open symbols correspond to He- 
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FIG 7.2. ^-extrapolation of several observables for 54 Fe calculated with the Kuo-Brown 
interaction KB3. 
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FIG 7.3. Comparison of SMMC results for selected observables with those obtained by 
conventional diagonalization |27J]. The calculations have been performed for 20 Ne within 
the sd model space. A pairing+quadrupole interaction, which is free from the sign problem 
discussed in section 7.2, has been used. 



(see subsection 7.2). 

Apart from demonstrating convergence in the various calculational parameters (A/3, num- 
ber of samples, etc.), SMMC algorithms have been validated by comparison with direct diag- 
onalization. Several examples of both the zero temperature and thermal SMMC formalisms 



are given in Refs. |26| , |2T|j . These calculations have been performed for Hamiltonians pos- 
sessing good sign properties (i.e. all V a < 0). For ground state and thermal properties, as 
well as for strength distributions, excellent agreement between direct diagonalization and 



the SMMC methods has been demonstrated p7| . Some examples are summarized in Fig. 
7.3. 

The ^-extrapolation has been tested for several nuclei in the s<i-shell |29| and in the low pf- 
shell [p5| . For example, using the parameter x — 4 SMMC can reproduce [55|] all observables 



(ground state energy, total B(E2), B(M1), GT strength etc.) for those nuclei in the mass 
range A = 48 — 56 for which direct diagonalization studies [51| or estimates on the basis of a 
series of truncation calculations [52~| exist. An example for the nuclei 54 ' 56 Fe is given in Fig. 
2.2. 

The extrapolation procedure in thermal SMMC studies has been validated for 44 Ti by 
comparison to a direct diagonalization calculation. The code CRUNCHER and the residual 
FPVH interaction were used to calculate the lowest 750 eigenstates up to an excitation 
energy of 43 MeV, which is sufficient to calculate the energy expectation value as a function 
of temperature up to T = 2 MeV. We note that exact reproduction of the CRUNCHER 
results was obtained only after a A/3 — > extrapolation which lowers the energies slightly 
compared to the calculation with finite A/3 = 1/32 MeV -1 . Fig. 7.4 compares the (A/?- and 
^-extrapolated) SMMC results with those obtained by direct diagonalization. Note that we 
cannot validate thermal properties requiring the calculation of transition matrix elements, 
as the latter are not numerically practical in direct diagonalization for thermal ensembles. 
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FIG 7.4. Comparison of SMMC results for the temperature dependence of the internal 
energy in 44 Ti with those obtained by conventional diagonalization code CRUNCHER. 



8 Selected results 

The previous sections have outlined the motivation, formalism, implementation, and valida- 
tion of Monte Carlo methods for treating the nuclear shell model. Of course, the ultimate 
goal of this work is to gain insight into the properties of real nuclei under a variety of con- 
ditions. In this section, we present a sampling of such results and their interpretation. The 
calculations presented should be viewed only as an indication of the power and potential 
of the method; undoubtedly, more work of this sort, abetted by ever- increasing computer 
power, will follow in the future. 



8.1 Ground state properties of j^f-shell nuclei 

While complete Ohu calculations can be carried out by direct diagonalization in the p- and 
s<i-shells, the exponentially increasing number of configurations limits such studies in the 



next (pf) shell to only the very lightest nuclei ]38], |51j. SMMC techniques allow calculation 
of groundstate observables in the full Ohu model space for nuclei throughout the pf -shell. 
Here, we discuss a set of such calculations that use the modified KB3 interaction jl^fl ; the 
single-particle basis is such that N s = 20 for both protons and neutrons. A more detailed 



description of the calculations and their results are found in Ref. . 

These studies were performed for 28 even-even Ti, Cr, Fe, Ni, and Zn isotopes; extension 
to heavier nuclei is dubious because of the neglect of the g 9 / 2 orbital. We have used (3 = 2 
MeV" 1 with fixed A/3 = 1/32 MeV" 1 ; test calculations at varying (3 show that the results 
accurately reflect the ground state properties. Some 4000-5000 independent Monte Carlo 
samples were taken for six equally-spaced values of g between —1 and 0, and observables 
were extrapolated to g — 1 using linear or quadratic functions, as required. We chose the 
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parameter x = 4, as described in Section 7.2. 

Figure 8.1 shows systematic results for the mass-defects, obtained directly from (H). The 
SMMC results have been corrected for the Coulomb energy, which is not included in the 



KB3 interaction, using 51 



H Coul = ^ U ■ 0.35 - ttz/0.05 + 7T • 7.289. (8.1) 

where ir and v are the numbers of valence protons and neutrons, respectively, and the energy 
is in MeV. As in Ref. [(jIH , we have increased the calculated energy expectation values by 
0.014-n(n— 1) MeV (where n = n+u is the number of valence nucleons) to correct for a "tiny" 
residual monopole defect in the KB3 interaction. In general, there is excellent agreement; 
the average error for the nuclei shown is +0.45 MeV, which agrees roughly with the internal 
excitation energy of a few hundred keV expected in our finite-temperature calculation. 

Figure 8.2 shows the calculated total E2 strengths for selected p/-shell nuclei. This quan- 
tity is defined as 

B(E2) = ((e p Q p + e n Q n ) 2 ) , (8.2) 

with 

Qp(n)=£ r i y 2(0i,&), (8-3) 
% 

where the sum runs over all valence protons (neutrons). The effective charges were chosen 
to be e p = 1.35 and e n = 0.35, while we used b = 1.01A 1 / 6 fm for the oscillator length. 
Shown for comparison are the B(E2) values for the 0^ — > 2\ transition in each nucleus; in 
even-even nuclei some 20-30% of the strength comes from higher transitions. The overall 
trend is well reproduced, For the nickel isotopes 58 > 60 > 64 Ni, the total B(E2) strength is known 
from (e, e') data and agrees very nicely with our SMMC results. 

The Gamow- Teller (GT) properties of nuclei in this region of the periodic table are crucial 
for supernova physics []TT| . These strengths are defined as 

B(GT ± ) = <(ar ± ) 2 >. (8.4) 

Note that the Ikeda sum rule, _B(GT_) — B(GT + ) = 3(N — Z) is satisfied exactly by our 
calculations. In Figure 8.3, we show the calculated B(GT + ) for all p/-shell nuclei for which 
this quantity has been measured by (n,p) reactions. (For the odd nuclei 51 V, 55 Mn and 
59 Co the SMMC calculations have been performed at (5 = 1 MeV -1 to avoid the odd- A sign 
problem. For even-even nuclei the GT strength calculated at this temperature is still a good 
approximation for the ground state value, as is shown below). Since the experimental values 
are normalized to measured /3-decay ft values, and it is believed that the axial coupling 
constant is renormalized from its free value = 1.26) to g& = 1 in nuclei |5"I| , we have 
multiplied all calculated results by (1/1. 26) 2 . Note that this results in generally excellent 
agreement between the calculations and the data. Thus, we support both the statement 
that gA is quenched to 1 in nuclei and the statement that complete shell-model calculations 
can account for the GT strength observed experimentally. Note that the quenching of g^ 
in the nuclear medium is not quite understood yet. It is believed to be related either to 



a second-order core polarization caused by the tensor force [^g] or to the screening of the 
Gamow- Teller operator by A-hole pairs [0 
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FIG 8.1. Upper panel (a): Comparison of the mass excesses AM as calculated within the 
SMMC approach with data. Lower panel (b): Discrepancy between the SMMC results for 
the mass excesses and the data, 5AM. The solid line shows the average discrepancy, 450 
keV, while the dashed lines show the rms variation about this value (from |55|). 
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FIG 8.2. Comparison of the experimental B(E2, Of — > 2+) strengths with the total B(E2) 
strength calculated in the SMMC approach for various p/-shell nuclei having either proton 
or neutron number of 28. For the nickel isotopes 58 > 60 > 64 Ni the total B(E2) strength (full 
squares) is known from inelastic electron scattering data (from |55| ). 
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FIG 8.4. Unrenormalized total Gamow- Teller strengths for various isotope chains as a func- 
tion of neutron holes (iV val — 20) in the p/-shell (from ]53|). 



The calculated systematics of the B(GT + ) strength for isotope chains throughout this 
region are shown in Figure 8.4, where it can be seen that GT + increases as the number 
of valence neutron holes increases (unblocking), and as the number of valence protons in- 
creases. These systematics were summarized phenomenologically from the experimental data 
in Ref. |j7| as 

B{GT+) = 0.045 -Z val - (20 - N^) , (8.5) 

which is derived by assuming that the residual interaction is strong enough to destroy all 
sub-shell structure, so that the entire p/-shell behaves as one large orbital. Figure 8.5 shows 
the linearity of B(GT + )/Z with (iV va i — 20) for the experimental data. 
The magnetic dipole strengths are defined as 

B(M1) = ((J2m{giUg s s}) 2 ). (8.6) 

i 

where /ijy is the nuclear magneton and gi,g s are the free gyromagnetic ratios for angular 
momentum and spin, respectively (gi = l,g s = 5.586 for protons, and g\ — 0, g s = —3.826 
for neutrons). A complete list of the total B(M1) strengths as calculated in SMMC studies 
of pf -shell nuclei is given in Ref. 



Brown and Wildenthal [21| suggested that, as a general rule, the strength of spin operators 



is renormalized in nuclei. To account for this renormalization, it is customary to use effective 



gyromagnetic ratios for the spin, g s = g s /1.26 in shell model calculations [21, 51]. If we follow 
this presciption, we calculate total B(M1) strengths for the N = 28 isotones (in units of 
H 2 N ) of 7.9 ± 0.7 ( 50 Ti), 11.7 ± 1.5 ( 52 Cr), and 10.2 ± 2 ( 54 Fe). For these nuclei, Richter 
and collaborators |7[J have used high-resolution electron scattering to study the B(M1) 
strength in an energy window large enough to contain most of the strength. They find 
B(M1) = 4.5 ± 0.5 ( 50 Ti), 8.1 ± 0.8 ( 52 Cr) and 6.6 ± 0.4 ( 54 Fe) for excitation energies 
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FIG 8.5. Experimental total Gamow- Teller strengths p2|-|55|, divided by the number of 
valence protons, as a function of neutron holes in the p/-shell. 



between 7 and 12 MeV [fOfl. If one considers that this energy window should contain about 
75% of the total strength [J7[J, our SMMC results appear to be consistent with the data. 



8.2 Pair structure of the nuclear ground state 

Numerous phenomenological descriptions of nuclear collective motion describe the nuclear 
ground state and its low-lying excitations in terms of bosons. In one such model, the Inter- 
acting Boson Model (IBM), L = (S) and L = 2 (D) bosons are identified with nucleon pairs 
having the same quantum numbers JHJ , and the ground state can be viewed as a condensate 
of such pairs. Shell model studies of the pair structure of the ground state and its variation 
with the number of valence nucleons can therefore shed light on the validity and microscopic 
foundations of these boson approaches. 

For many purposes, it appears sufficient to study the BCS pair structure in the ground 
state. The BCS pair operator for protons can be defined as 



jm>0 



mPjm 5 



17) 



where the sum is over all orbitals with m > and = {—) J ""Pj m 



'•' rn ^- is the time-reversed 
operator. Thus, the observable A^ A and its analog for neutrons are measures of the numbers 
of J = 0,T = 1 pairs in the groundstate. For an uncorrelated Fermi gas, we have simply 



A) = £ 



11; 



2(2.7 + 1) 



(8.8) 
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where the rij = (p\ rn Pj m ) are the occupation numbers, so that any excess of (At A) in our 
SMMC calculations over the Fermi-gas value indicates pairing correlations in the ground- 
state. 

Fig. 8.6 shows the SMMC expectation values of the proton and neutron BCS-like pairs, 
obtained after subtraction of the Fermi gas value (8.8), for three chains of isotopes. As 
expected, these excess pair correlations are quite strong and reflect the well-known coherence 
in the ground states of even-even nuclei. Note that the proton BCS-like pairing fields are 
not constant within an isotope chain, showing that there are important proton-neutron 
correlations present in the ground state. The shell closure at N = 28 is manifest in the 
neutron BCS-like pairing. As is demonstrated in Fig. 8.7, the proton and neutron occupation 
numbers show a much smoother behavior with increasing A. 

It should be noted that the BCS form (8.7) in which all time-reversed pairs have equal 
amplitude is not necessarily the optimal one and allows only the study of S-pair structure. 
To explore the pair content of the ground state in a more general way, we define proton pair 
creation operators 

^j,(jaJ b ) = / J_ [at | x a) b ]j, . (8.9) 

These operators are boson-like in the sense that 

[A^(j a j b ),Aj,(j a j b )} = 1 + 0(n/2j + 1) ; (8.10) 

i.e., they satisfy the expected commutation relations in the limit of an empty shell. 
We construct bosons JfX as 

BL, = E i>ax(jaj b )A[^j a j b ) , (8.11) 

jajb 

where a — 1, 2, . . . labels the particular boson and the "wave function" ip satisfies 

2 ^*ajUajb)^f3j(jajb) = 5 a p • (8.12) 

jajb 

(Note that ip is independent of /x by rotational invariance.) 

To find ip and n aJ = I^X-Btj _B aJ/1 ), the number of bosons of type a and multipolarity J, 
we compute the quantity H^{^j^{jaib)^-jfA,{jcjd)) , which can be thought of as an hermitian 
matrix M^ a , in the space of orbital pairs (j a jb)] its non- negative eigenvalues define the n aJ 
(we order them so that n\j > n 2 j > ■ ■ .), while the normalized eigenvectors are the ip a j(jajb)- 
The index a distinguishes the various possible bosons. For example, in the complete pj-shell 
the square matrix M has dimension Nj = 4 for J = 0, Nj = 10 for J = 1, Nj = 13 for 
J = 2,3. 

The presence of a pair condensate in a correlated ground state will be signaled by the 
largest eigenvalue for a given J, riu, being much greater than any of the others; tpu will 
then be the condensate wavefunction. In Fig. 8.8 we show the pair matrix eigenvalues n a j 
for the three isovector J = + pairing channels as calculated for the iron isotopes 54_58 Fe. We 
compare the SMMC results with those of an uncorrelated Fermi gas, where we can compute 
(A^A) using the factorization 

(ala^a^as) = npn a (6fa6as - 5 ps 5 ai ) , (8.13) 
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FIG 8.6. SMMC expectation values of proton and neutron BCS-like pairs after subtraction 
of the Fermi gas value (from f55[). 
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where the rip = (alap) are the occupation numbers. Additionally, Fig. 8.8 shows the 
diagonal matrix elements of the pair matrix M aa . As expected, the protons occupy mainly 
fj/2 orbitals in these nuclei. Correspondingly, the (A^A) expectation value is large for this 
orbital and small otherwise. For neutrons, the pair matrix is also largest for the /V/ 2 orbital. 
The excess neutrons in 56 > 58 Fe occupy the p 3 / 2 orbital, signalled by a strong increase of the 
corresponding pair matrix element M 2 2 in comparison to its value for 54 Fe. Upon closer 
inspection we find that the proton pair matrix elements are not constant within the isotope 
chain. This behavior is mainly caused by the isoscalar proton-neutron pairing. The dominant 
role is played by the isoscalar 1 + channel, which couples protons and neutrons in the same 
orbitals and in spin-orbit partners. As a consequence we find that, for 54 ' 56 Fe, the proton 
pair matrix in the f§/2 orbital, M 33 , is larger than in the p^/2 orbital, although the latter is 
favored in energy. For 58 Fe, this ordering is inverted, as the increasing number of neutrons 
in the p%/2 orbital increases the proton pairing in that orbital. 

After diagonalization of M, the J = proton pairing strength is essentially found in one 
large eigenvalue. Furthermore we observe that this eigenvalue is significantly larger than the 
largest eigenvalue on the mean- field level (Fermi gas), supporting the existence of a proton 
pair condensate in the ground state of these nuclei. The situation is somewhat different for 
neutrons. For 54 Fe only little additional coherence is found beyond the mean-field value, 
reflecting the closed-subshell neutron structure. For the two other isotopes, the neutron 
pairing exhibits two large eigenvalues. Although the larger one exceeds the mean-field value 
and signals noticeable additional coherence across the subshells, the existence of a second 
coherent eigenvalue shows the shortcomings of the BCS-like pairing picture. 

It has long being anticipated that J = + proton-neutron correlations play an important 
role in the ground states of N = Z nuclei. To explore these correlations we have performed 
SMMC calculations of the N = Z nuclei in the mass region A = 48 — 56. Note that for these 
nuclei the pair matrix in all three isovector + channels essentially exhibits only one large 
eigenvalue, related to the /V/2 orbital. We will use this eigenvalue as a convenient measure 
of the pairing strength. As the even-even N = Z nuclei have isospin T = 0, the expectation 
values of A*A are identical in all three isovector + pairing channels. This symmetry does not 
hold for the odd-odd N = Z nuclei in this mass range, which have T = 1 ground states, and 
(A' A) can be different for proton-neutron pairs than for like-nucleons pair (the expectation 
values for proton pairs and neutron pairs are identical). We find the proton-neutron pairing 
strength significantly larger for odd-odd N = Z nuclei than in even-even nuclei, while the + 
proton and neutron pairing shows the opposite behavior, in both cases leading to an odd-even 
staggering, as displayed in Fig. 8.9. This staggering is caused by a constructive interference 
of the isotensor and isoscalar parts of A^A in the odd-odd N = Z nuclei, while they interfere 
destructively in the even-even nuclei. The isoscalar part is related to the pairing energy, and 
is found to be about constant for the nuclei studied here. 



8.3 Thermal properties of pf-sheU nuclei 

The properties of nuclei at finite temperatures are of considerable experimental W% [72] and 



theoretical interest []73], p4\ . They are clearly quite important in various astrophysical scenar- 
ios. For example, electron capture on nuclei plays an essential role in the early presupernova 
collapse |TT|. In that context, it is important to know the temperature dependence of the 
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FIG 8.8. Content of isovector + pairs and isoscalar 1 + pairs in the ground states of the 
isotopes 54_58 Fe. The upper panel shows the diagonal matrix elements of the pair matrix 
M aa . The index a = 1, .., 4 refers to + pairs in the ff/ 2 , Villi fs/2, and pi/ 2 orbitals, 
respectively. For the isoscalar pairs a = 1,2, 3 refers to (/V/2) 2 , (jV/2/5/2) and (fs/2) 2 pairs, 
respectively. The middle panel gives the eigenvalues of the pair matrix; for the isoscalar 
pairs only the 3 largest are shown. The lower panel gives the eigenvalues of the pair matrix 
for the uncorrelated Fermi gas case using Eq. (8.13) (from J77[ ) . 
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Gamow- Teller strength. 

Theoretical studies of nuclei at finite temperature have been based mainly on mean-field 
approaches, and thus only consider the temperature dependence of the most probable con- 
figuration in a given system. These approaches have been criticized due to their neglect 
of quantum and statistical fluctuations f75| . The SMMC method does not suffer this de- 
fect and allows the consideration of model spaces large enough to account for the relevant 
nucleon-nucleon correlations at low and moderate temperatures. 

We have performed SMMC calculations of the thermal properties of several even-even 
nuclei in the mass region A = 50 — 60 [76, 77j. As a typical example, we discuss in the 



following our SMMC results for the nucleus 54 Fe, which is very abundant in the presupernova 
core of a massive star. 

Our calculations include the complete set of 1^3/2,1/20/7/2,5/2 states interacting through the 
realistic Brown- Richter Hamiltonian [|3£|. (SMMC calculations using the modified KB3 inter- 
action give essentially the same results.) Some 5 x 10 9 configurations of the 8 valence neutrons 
and 6 valence protons moving in these 20 orbitals are involved in the canonical ensemble. 
The results presented below have been obtained with a time step of A/3 = 1/32 MeV -1 
using 5000-9000 independent Monte Carlo samples at seven values of the coupling constant 
g spaced between —1 and and the value x — 4. 

The calculated temperature dependence of various observables is shown in Fig. 8.10. In 
accord with general thermodynamic principles, the internal energy U steadily increases with 
increasing temperature ||76|| . It shows an inflection point around T « 1.1 MeV, leading to 
a peak in the heat capacity, C = dll/dT, whose physical origin we will discuss below. The 
decrease in C for T > 1.4 MeV is due to our finite model space (the Schottky effect @); 
we estimate that limitation of the model space to only the p/-shell renders our calculations 
of 54 Fe quantitatively unreliable for temperatures above this value (internal energies U 
15 MeV). The same behavior is apparent in the level density parameter, a = C/2T. The 
empirical value for a is A/8 MeV = 6.8 MeV" 1 which is in good agreement with our results 
for T ~ 1.1—1.5 MeV. 

We also show in Fig. 8.10 the expectation values of the BCS-like proton-proton and 
neutron-neutron pairing fields, (At A). At low temperatures, the pairing fields are signif- 
icantly larger than those calculated for a non-interacting Fermi gas, indicating a strong 
coherence in the ground state. With increasing temperature, the pairing fields decrease, and 
both approach the Fermi gas values for T « 1.5 MeV and follow it closely for even higher 
temperatures. Associated with the breaking of pairs is a dramatic increase in the moment of 
inertia, / = (J 2 )/3T, for T = 1.0-1.5 MeV; this is analogous to the rapid increase in mag- 
netic susceptibility in a superconductor. At temperatures above 1.5 MeV, / is in agreement 
with the rigid rotor value, 10.7^ /MeV; at even higher temperatures it decreases linearly 
due to our finite model space. 

In Fig. 8.11, we show various static observables. The Ml strength unquenches rapidly 
with heating near the transition temperature. However, for T = 1.3-2 MeV B(M1) re- 
mains significantly lower than the single-particle estimate (41 suggesting a persistent 
quenching at temperatures above the like-nucleon depairing. This finding is supported by the 
near-constancy of the Gamow- Teller f3 + strength, B(GT + ), for temperatures up to 2 MeV. 
As the results of Ref. |77j demonstrate that neutron-proton correlations are responsible for 
much of the GT quenching in iron nuclei at zero temperature, we interpret the present re- 
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suits as evidence that isovector proton-neutron correlations persist to higher temperatures 
(see below). We have verified that, in our restricted model space, both the GT + and Ml 
strengths unquench at temperatures above 2 MeV and that, in the high-temperature limit, 
they both approach the appropriate Fermi gas values. We note that it is often assumed 
in astrophysical calculations that the GT strength is independent of temperature our 
calculations demonstrate that this is true for the relevant temperature regime (T < 2 MeV). 
We also note that a detailed examination of the occupation numbers of the various orbitals 
show no unusual variation as the pairing vanishes. 

Although the results discussed above are typical for even-even nuclei in this mass region 
(including the N = Z nucleus 52 Fe), they are not for odd-odd N = Z nuclei. This is 
illustrated in Fig. 8.12 which shows the thermal behavior of several observables for 50 Mn 
(N = Z = 25), calculated in a SMMC study within the complete p/-shell using the KB3 
interaction. In contrast to even-even nuclei, the total Gamow- Teller strength is not constant 
at low temperatures, but increases by about 50% between T = 0.4 MeV and 1 MeV. The 
5 (Ml) strength decreases significantly in the same temperature interval, while for even- 
even nuclei, it increases steadily. A closer inspection of the isovector J = and isoscalar 
J = 1 pairing correlations holds the key to the understanding of these differences. The J = 
isovector correlations are studied using the BCS pair operators (8.6), with a similar definition 
for proton-neutron pairing. For the isoscalar J = 1 correlations we have interpreted the trace 
of the pair matrix M J=1 as an overall measure for the pairing strength, 

Psm = EH = Y, M L (8-14) 

Note that at the level of the non-interacting Fermi gas, proton-proton, neutron-neutron 
and proton-neutron J = correlations are identical for N = Z nuclei. However, the residual 
interaction breaks the symmetry between like-pair correlations and proton-neutron corre- 
lations in odd-odd N = Z nuclei. As is obvious from Fig. 8.12, at low temperatures 
proton-neutron pairing dominates in 50 Mn, while pairing among like nucleons shows only a 
small excess over the Fermi gas values, in strong contrast to even-even nuclei. 

A striking feature of Fig. 8.12 is that the isovector proton-neutron correlations decrease 
strongly with temperature and have essentially vanished at T = 1 MeV, while the isoscalar 
pairing strength remains about constant in this temperature region (as it does in even-even 
nuclei) and greatly exceeds the Fermi gas values. We also note that the pairing between like 
nucleons is roughly constant at T < 1 MeV. The change of importance between isovector 
and isoscalar proton-neutron correlations with temperature is nicely reflected in the isospin 
expectation value, which decreases from < T 2 >= 2 at temperatures around 0.5 MeV, 
corresponding to the dominance of isovector correlations, to < T 2 >= at temperature 
T = 1 MeV, when isoscalar proton- neutron correlations are most important. The low- 
temperature behavior of the Gamow- Teller and B(M1) strength is related to the fading of 
the isoscalar proton-neutron correlations. For example, the Gamow- Teller strength increases 
as transitions between the same orbital (mainly /V/2) are less quenched by weaker isovector 
proton-neutron correlations. 
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FIG 8.10. Temperature dependence of various observables in 54 Fe. Monte Carlo points with 
statistical errors are shown at each temperature T. In the left-hand column, the internal 
energy, U, is calculated as (H) — E , where H is the many-body Hamiltonian and E the 
ground state energy. The heat capacity C is calculated by a finite-difference approximation 
to dU/dT, after U(T) has been subjected to a three-point smoothing, and the level density 
parameter is a = C/2T. In the right-hand column, we show the expectation values of the 
squares of the proton and neutron BCS pairing fields. For comparison, the pairing fields 
calculated in an uncorrelated Fermi gas are shown by the solid curve. The moment of inertia 
is obtained from the expectation values of the square of the total angular momentum by 
/ = (3(j 2 )/3 (from |76|). 
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FIG 8.11. The upper left panel shows, for 54 Fe, the total magnetic dipole strength, B(M1) in 
units of nuclear magnetons; it is calculated using free-nucleon ^-factors. The lower left panel 
shows the GT (3 + strength, while the upper and lower right panels show the isoscalar (Q = 
Q P + Qn) and isovector {Q v = Q p — Q n ) quadrupole strengths. In these latter observables, 
the quadrupole operators are r 2 Y 2 , and results are given in terms of the oscillator length, 
b = 1.96 fm (from |7F 
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FIG 8.12 Temperature dependence of various observables in 50 Mn. The left panels show 
(from top to bottom) the total energy, the moment of inertia, the B(M1) strength, and the 
Gamow- Teller strength, while the right panels exhibit the expectation values of the isospin 
operator (T 2 ), the isovector J = proton-proton and proton-neutron BCS pairing fields, 
and the isoscalar J = 1 pairing strength, as defined in the text. For comparison, the solid 
lines indicate the Fermi gas values of the BCS pairing fields and J = 1 pairing strength. 
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8.4 Pair correlations and thermal response 

All SMMC calculations of even-even nuclei in the mass region A = 50 — 60 show that the 
BCS-like pairs break at temperatures around 1 MeV. Three observables exhibit a particularly 
interesting behavior at this phase transition: a) the moment of inertia rises sharply; b) 
the Ml strength shows a sharp rise, but unquenches only partially; and c) the Gamow- 
Teller strength remains roughly constant (and strongly quenched). Note that the B(M1) 
and B(GT + ) strengths unquench at temperatures larger than rs 2.6 MeV and in the high- 
temperature limit approach the appropriate values for our adopted model space. 
Ref. [77j has studied the pair correlations in the four nuclei 54_58 Fe and 56 Cr for the various 



isovector and isoscalar pairs up to J = 4, tentatively interpreting the sum of the eigenvalues 
of the matrix M J (8.14) as an overall measure for the pairing strength. Note that the pairing 
strength, as defined in (8.12), is non-zero at the mean-field level. The physically relevant 
pair correlations P^ orr are then defined as the difference of the SMMC and mean-field pairing 
strengths. 

Detailed calculations of the pair correlations have been performed for selected temperatures 
in the region between T = 0.5 MeV and 8 MeV. Fig. 8.13 shows the temperature dependence 
of those pair correlations that play an important role in understanding the thermal behavior 
of the moment of inertia and the total Ml and Gamow- Teller strengths. 

The most interesting behavior is found in the J = proton and neutron pairs. There is 
a large excess of this pairing at low temperatures, reflecting the ground state coherence of 
even-even nuclei. With increasing temperature, this excess diminishes and vanishes at around 
T = 1.2 MeV. We observe further from Fig. 8.13 that the temperature dependence of the 
J = proton-pair correlations are remarkably independent of the nucleus, while the neutron 
pair correlations show interesting differences. At first, the correlation excess is smaller in the 
semimagic nucleus 54 Fe than in the others. When comparing the iron isotopes, the vanishing 
of the neutron J = correlations occurs at higher temperatures with increasing neutron 
number. 

The vanishing of the J = proton and neutron pair correlations is accompanied by an 
increase in the correlations of the other pairs. For example, the isovector J = proton- 
neutron correlations increase by about a factor 3 after the J = proton and neutron pairs 
have vanished. The correlation peak is reached at higher temperatures with increasing 
neutron number, while the peak height decreases with neutron excess. 

The isoscalar proton-neutron J = 1 pairs show an interesting temperature dependence. At 
low temperatures, when the nucleus is still dominated by the J = proton and neutron pairs, 
the isoscalar proton-neutron correlations show a noticeable excess but, more interestingly, 
they are roughly constant and do not directly reflect the vanishing of the J = proton and 
neutron pairs. However, at T > 1 MeV, where the J = proton- and neutron-pairs have 
broken, the isoscalar J = 1 pair correlations significantly increase and have their maximum 
at around 2 MeV, with peak values of about twice the correlation excess in the ground state. 
In contrast to the isovector J = proton-neutron pairs, the correlation peaks occur at lower 
temperatures with increasing neutron excess. We also observe that these correlations fade 
rather slowly with increasing temperature. 

We note that the SMMC results for the moment of inertia and the B(M1) and B(GT + ) 
strengths are significantly smaller than the mean-field values. The quenching of these ob- 
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FIG 8.13 Pair correlations, as defined in Eq. (8.13), for isovector + and isoscalar 1 + pairs 
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FIG 8.14. The temperature dependence of the quenching of the moment of inertia and of 
the total B(M1) and B{GT + ) strengths (see text) is compared to the one of selected pairing 
correlations (in arbitrary units). The results are shown for 56 Fe (from ||77|| ). 
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servables, defined as the difference of the mean-field and SMMC results, can be traced back 
to nucleon correlations beyond the mean-field. This is demonstrated in Fig. 8.14, exempli- 
fied for 56 Fe. From its definition (J ~ (J 2 ) = (( J n + J p ) 2 ), where J p , J n are the proton and 
neutron angular momenta, respectively), one expects that the moment of inertia is sensitive 
to the proton, neutron and proton- neutron correlations. In fact, the sum of the isovector 
J = pairing correlations shows nearly the same temperature behavior as the moment of 
inertia quenching. The drop in the quenching at T « 1 MeV clearly signals the pairing 
phase transition observed in the J = proton and neutron correlations (see Fig. 8.13). The 
temperature dependence of the Ml quenching is well described by that of the sum of the 
J = proton and isoscalar J = 1 proton-neutron correlations. These two correlations reflect 
the two components in the Ml strength. The orbital part is sensitive to the J = proton 
pairing correlations (the gyromagnetic moment g\ is zero for neutrons), while the quenching 
of the spin component is dominated by isoscalar proton-neutron correlations. 

The quenching of the B(GT + ) strength is roughly constant for T < 1 MeV. It then in- 
creases slightly, before it slowly dies out following a maximum at T m 2.5 MeV. While the 
temperature dependence of the Gamow- Teller strength is driven by the isoscalar proton- 
neutron correlations, the overall measure of the pairing correlations, as applied here, is too 
simple for a quantitative reproduction. 



8.5 Gamow- Teller quenching in the 100 Sn region 



The double-magic nucleus 
Ganil HBO, RT1 



100 Sn has recently been observed for the first time at GSI and 
Detailed studies of this nucleus promise further insight into the shell structure 
in proton-rich nuclei and in particular into nature of Gamow- Teller (GT) quenching. The 
reasoning is as follows: the shell gap between the g 9 / 2 -orbital and the other orbitals in the 
0g-ld-2s shell is expected to support a closed-shell structure for both neutrons and protons 
and suppress correlations involving the higher orbitals. The valence protons, occupying g 9 / 2 - 
orbitals, will then be transformed with relatively less hindrance into (77/2-neutrons by the ar + 
operator. Thus the quenching of the GT + -strength due to configuration mixing is expected 
to be small. 

Despite this reasoning, the observed quenching of the GT strength in the N = 50 isotones 
close to 100 Sn like 94 Ru, 96 Pd and 98 Cd is surprisingly large Strongly truncated shell model 
calculations are unable to reproduce the strong quenching, as the configuration mixing in 
these restricted model spaces was found to be insufficient [S3L EM. The agreement between 
theory and experiment has been improved by renormalization of the GT operator to consider 
(first-order) corrections to the wave functions |[5|, |86|, 83|, jB4 |. As shell model studies of lighter 
nuclei clearly stressed the need of a complete Ohu shell model basis, it appears unlikely that 
first-order corrections are sufficient to recover the full quenching of the GT strength due to 
configuration mixing. 

Shell model Monte Carlo calculations have been performed for 94 Ru, 96 Pd, 96 ' 98 Cd, and 
100 Sn in the complete 0g-ld-2s major oscillator shell using a G-matrix two-body interaction 



The Coulomb-corrected SMMC mass 



built from the Paris nucleon-nucleon interaction 
excesses agree well with the latest mass compilation of Myers and Swiatecki |8q| . Table 3 
gives the SMMC results for the Gamow- Teller strength, renormalizing the GT operator by 
1/1.26. The SMMC calculation predicts a noticeable quenching of the GT strength in the 
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Table 3: Comparison of the B(GT) strengths and quenching factors q, as calculated in our 
SMMC approach, with the data. The single particle estimates (SPE) for the B(GT) strength 
are also given (from |87j| ). 



Nucleus 


B(GT+) (SPE) 


B(GT+) (SMMC) 


q (SMMC) 


q (expt) |82 


94 Ru 


7.11 


1.48±0.1 


4.96±0.3 


7.2±0.6 


96p d 


10.67 


2.68±0.1 


3.98±0.2 


4 fi +1 - 7 

4 -°~0.8 


96 Cd 


16.18 


5.40±0.1 


2.98±0.1 


? 


98 Cd 


14.22 


4.42±0.1 


3.22±0.1 


4 1 +1 


ioog n 


17.78 


6.50±0.1 


2.74±0.1 


? 



= 50 isotones, in accord with experiment. A detailed comparison with experiment is 
hampered by the fact that the measurement of the /3 + -decay can only determine the GT 
strength below a certain energy threshold, so that the experiment actually provides only 
upper bounds for the quenching. In the case of 96 Pd and 98 Cd it is expected that the 

while for 94 Ru 



measurements miss less than about 10% of the total GT strength |8^ 
a significant amount of GT strength is predicted at energies that are not accessible in f3 + 
experiments [^, |84|. It thus seems that SMMC calculations are in excellent agreement with 
the observed GT quenching for the nuclei 96 Pd and 98 Cd. For 94 Ru the SMMC result appears 
to be in qualitative agreement with observation. However, a quantitative comparison requires 
a detailed knowledge of the GT strength distribution in the daughter nucleus 94 Tc, which is 
outside the scope of the calculation of Ref. . 



The SMMC calculations |87j also predict a noticeable quenching of the GT strength in the 
double-magic nucleus 100 Sn. Here, the total B(GT + )=6.5±0.1 should be compared to the 
single particle estimate of 17.78, corresponding to a quenching factor q = 2.74 ± 0.1. This 
quenching factor is larger than predicted by recent truncated shell model calculations 
However, the ratio of quenching factors for 98 Cd and 100 Sn is found to be less interaction de- 
pendent than the actual values |B1. The SMMC calculation yields q( 98 Cd)/q( 100 Sn)=1.18± 



0.07, within the range 1.23 - 1.3 given in [R4 



8.6 Electron capture and presupernova collapse 

The core of a massive star at the end of hydrostatic burning is stabilized by electron degener- 
acy pressure as long as its mass does not exceed the appropriate Chandrasekhar mass M C h- 



If the core mass exceeds Mqh, electrons are captured by nuclei ||TTJ to avoid a violation of 
the Pauli principle: 

e~ + (Z,A) (Z- I, A) + v e (8.15) 

The neutrinos can still escape from the core, carrying away energy. This is accompanied by 
a loss of pressure and the collapse is accelerated. 
For many of the nuclei that determine the electron capture rate in the early stage of the 



presupernova |79[, Gamow- Teller (GT) transitions contribute significantly to the electron 
capture rate. Due to insufficient experimental information, the GT + transition rates have so 
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far been treated only qualitatively in presupernova collapse simulations, assuming the GT + 
strength to reside in a single resonance, whose energy relative to the daughter ground state 



has been parametrized phenomenologically Rl; the total GT + strength has been taken from 



the single particle model. Recent (n,p) experiments p^]-[p^|, however, show that the GT + 
strength is fragmented over many states, while the total strength is significantly quenched 
compared to the single particle model (see section 8.1). (A recent update of the GT + rates 
for use in supernova simulations assumed a constant quenching factor of 2 JFjJ ) . 

In a series of truncated shell model calculations, Aufderheide and collaborators have 
demonstrated that a strong phase space dependence makes the Gamow- Teller contributions 
to the presupernova electron capture rates more sensitive to the strength distribution in the 



daughter nucleus than to the total strength In this work it also became apparent that 



complete Ohu studies of the GT + strength distribution are desirable. Such studies are now 
possible using the SMMC approach. 

To determine the GT + strength distribution we have calculated the response function of 
the o~t + operator, R GT (r), as defined in Eq. (3.6). As the strength function Sqt(E) is 
the inverse Laplace transform of Rgt{t), we have used the Maximum Entropy technique, 
described in section 6.4, to extract Sqt(E). For the default model m(E) in (6.21) we adopted 
a Gaussian whose centroid is given by the slope of IoRgt(t) at small r and whose width is 
2 MeV. 

As first examples we have studied several nuclei ( 51 V, 54 ' 56 Fe, 58 ' 60 ' 62 Ni, and 59 Co), for 
which the Gamow- Teller strength distribution in the daughter nucleus is known from (n,p) 
experiments |]62|]-||66|]. Note that the electron capture by these nuclei, however, plays only a 



minor role in the presupernova collapse. As SMMC calculates the strength function within 
the parent nucleus, the results have been shifted using the experimental Q-value. The 
Coulomb correction has been performed using Eq. (8.1). For all nuclei, the SMMC approach 
calculates the centroid and width of the strength distribution in good agreement with data. 
(Following [52] the calculated strength distributions have been folded with Gaussians of 
width 1.77 MeV to account for the experimental resolution.) Fig. 8.15 shows the response 
function Rgt{t) for 54 Fe and 59 Co and compares the extracted GT + strength distribution 
with data. (As discussed in sections 8.1 and 8.5, the Gamow- Teller operator has been 
renormalized by the factor 0.77). The centroid of the GT + strength distributions is found to 
be nearly independent of temperature (calculations for odd nuclei could be performed only 
in the temperature range T = 0.8 — 1.2 MeV), while its width increases with temperature. 

Following the formalism described in Refs. |89|. [n| , the Gamow- Teller contributions to the 
electron capture rates under typical presupernova conditions have been calculated assuming 
that the electrons obey a Fermi-Dirac distribution with a chemical potential adopted from the 
stellar trajectory at the electron-to-nucleon ratio corresponding to the respective nucleus |79 |. 
The calculation have been performed for both the SMMC and experimental GT + strength 
distributions |6^j-[Q. The electron capture rates so obtained agree within a factor of two 
for temperatures T = (3 — 5) x 10 9 Kelvin, which is the relevant temperature regime in the 
presupernova collapse J7S[ (Fig. 8.16). Note that this level of agreement is as good as that 
obtained by Aufderheide et al. These authors calculate the GT + strength distribution 
within strongly truncated shell model studies by adjusting the single particle energies to 
fit the observed strength distribution and then normalize the total calculated strength to 
the measured B(GT + ) value As this approach obviously requires a knowledge of the 
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FIG 8.15. SMMC GT + response functions (left side) and GT + strength distributions (right 
side) for 54 Fe (upper panel), 59 Co (middle panel), and 55 Co (lower panel). The energies refer 
to the daughter nuclei. The dashed histograms show the experimental strength distribution 
as extracted from (n,p) data p3J. 



experimental strength distribution, the SMMC calculations in the complete pf shell are 
capable of predicting the desired strength distribution with a similar accuracy. Thus, it is 
for the first time possible to calculate with a reasonable accuracy the electron capture rate 
for those nuclei like 55 Co or 56 Ni which dominate the electron capture process in the early 
presupernova collapse ]7j|. The SMMC calculation places the GT + strength for 55 Co at a 
higher energy in the daughter nucleus than the phenemenological parametrization of Fuller 
et al. [jSOfl - Correspondingly, the electron capture rate on this nucleus, which is estimated to 
be dominated by the GT + contribution, is smaller than previously assumed. Studies of the 
effect of this smaller rate on the dynamics of the early presupernova collapse must await the 
completion of SMMC calculations for other important p/-shell nuclei now in progress. 



8.7 Temperature dependence of the nuclear symmetry energy 

The size of the core mantle (the difference between the iron and homologous cores) is a major 



determinant of the supernova mechanism [11]. It has been argued recently that the mantle 



size might be significantly smaller than generally calculated in supernova models because of 
an overlooked temperature dependence of the electron capture process [|91|]. The detailed 
reasoning is as follows ||92|| : When the nucleus is described by a static mean field model, 
dynamical effects associated with collective surface vibrations are conventionally embodied 
in an effective nucleon mass, m* (as are spatial non-localities in the mean field). Donati 
P2| studied the coupling of the mean field single-particle levels to the collective 



et al. 
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FIG 8.16. Gamow- Teller contributions to the electron capture rates on 54 Fe and 59 Co as a 
function of temperature (Tg measures the temperature in 10 9 Kelvin.) The calculations have 
been performed for the experimental strength distribution [63] and the SMMC result. 
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surface vibrations within the quasiparticle random phase approximation and found that m* 
decreases with increasing temperature for T < 2 MeV, the temperatures relevant for the 
presupernova collapse. They attributed this behavior to a reduction of collectivity at low 
excitation energies. In the Fermi gas model, this decrease of m* induces an increase in the 
symmetry energy contribution to the nuclear binding energy 



E sym (T) = 6 sym (T) ( N ^ , (8.16) 



where N, Z, and A are the neutron, proton, and mass numbers of the nucleus, respectively. 
Quantitatively, Donati et al. estimate that b sym (T) increases by about 2.5 MeV as T increases 
from to 1 MeV (6 sym (0) « 28 MeV §). Importantly, a larger value of E sym (T) would 
hinder electron capture in the presupernova environment, and thus reduce the size of the 
core mantle, so that the shock wave would need significantly less energy to stop the collapse 
and explode the star. 

To explore the temperature dependence of the symmetry energy, SMMC calculations were 
used |)3] to study the thermal properties of various pairs of isobars in the mass region 
A = 54 — 64 which is important for the presupernova collapse; this includes two of the three 
nuclei studied in Ref. @ ( 64 Zn and 64 Ni). 



The discussion focuses on the temperature dependence of the internal energy of a nucleus 
(with proton and neutron numbers Z and N), which is defined as 

U(N, Z,T) = (H)(T) - (H)(T = 0) , (8.17) 

where (H)(T) is the expectation value of the Hamiltonian in the canonical ensemble at 
temperature T. Guided by the semi-empirical parametrization of the binding energies (e.g., 
the Bethe-Weizsacker formula), the difference of U(T) for two isobars 

AU(T) = U(Ni,Zi,T)-U(N 2 ,Z 2 ,T) -N x + Z x = N 2 + Z 2 = A (8.18) 

should reflect the proposed temperature dependence of the symmetry energy contribution 
to the binding energies. This proportionality is valid only if the other contributions to the 
binding energy, particularly nuclear structure effects, are negligible in AU(T). To investigate 
this point, Dean et al. P3| performed calculations for several pairs of isobars for mass 



numbers A = 54 ( 54 Fe, 54 Cr), A = 56 ( 56 Fe, 56 Cr), A = 58 ( 58 Ni, 58 Fe), A = 60 ( 60 Ni, 60 Fe), 
A = 62 ( 62 Zn, 62 Ni) and A = 64 ( 64 Zn, 64 Ni). 

The experimental knowledge of the excitation spectra for the nuclei studied appeared to be 
sufficient to derive U(T) for T < 0.5 MeV solely from data. For higher temperatures these 
values were supplemented with the appropriate SMMC results for U(T) — U(T = 0.5 MeV). 
Fig. 8.17 compares AU(T) for T = 1 MeV as calculated by the combination of data and 



SMMC results [^3J with the predictions of Ref. . For all isobaric pairs with the exception 
of ( 64 Ni, 64 Zn), the values for AU(T) are compatible with zero (no increase of fe sy mm with 
temperature). However, we note that, in agreement with Ref. |92[|, the SMMC calculation 



shows an increase in 6 symm for the A = 64 pair ( 64 Ni, 64 Zn). However, it does not confirm 
that this increase is generic. 
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FIG 8.17. Comparison of the values for AU(T) predicted by Donati et al. (open circles, 
[P^l ) with the results calculated from experimental data and SMMC calculations for selected 
pairs of isobars, as described in the text (full circles). 



8.8 170 Dy at finite rotation and temperature 



A first SMMC calculation p5| to describe rare-earth nuclei used the Z = 50-82 shell for 
protons (2slc£0<77/ 2 0/iii/2) and the iV = 82-126 shell for neutrons (2plf0g 9 / 2 0ii 3 /2). The 
Hamiltonian chosen was of the pairing plus quadrupole form given by 



X 



Q-Q 



(8.19) 



where e a are the single particle energies, Ppi n yPp{n) are the monopole pair creation and 

annihilation operators for protons and neutrons, and Q = Q p +Q n is the quadrupole operator. 
The pairing strengths g p ( n ), the quadrupole interaction strength x, and the single particle 
energies used were taken from |40| . 



The nucleus considered was 170 Dy. It involved 16 valence protons and 22 valence neutrons 
moving among the N s = 32 and 44 single-particle states, respectively, giving a total of 10 21 
m-scheme determinants. The calculations used A/3 = 0.0625 MeV -1 and N t — 8 to 64 time 
slices. Due to limited computer memory and cycles at that time, it was necessary to perform 
canonical analyses of the fields generated through grand-canonical sampling; all other work 
in this Report uses canonical sampling. In particular, for the 170 Dy calculation, canonical 
observables (subscript "c") are given in terms of the grand-canonical sampling (subscript 
"G" ) as 



<n>« 



$DaW Ga <S> ca [Cca/CGa] 



(8.20) 



where ( a = TrU a . Although this is an exact expression, the fluctuations in [Cca/Cca] deter- 
mine how precise the Monte Carlo evaluation can be. We found that the fluctuations are 
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less than 10% (as the number dispersion in the grand-canonical ensemble is small), so that 
canonical observables can be calculated with good precision. 

In Fig. 8.18 we show the static observables (for the uncranked system) in both the grand- 
canonical and canonical formalisms. We calculated observables canonically, using grand 
canonical fields, up to (3 — 2.0 in order to demonstrate that, for these nuclear systems, either 
method may be used in this kind of calculation. As the temperature decreases, the nucleus 
becomes more deformed. Relaxation of the expectation values of H and J 2 is also clearly 
seen. The sign in these cases is identically one for all temperatures. 

Cranking calculations in which H — > H — ujJ z were also performed. The systematics 
are shown in Fig. 8.19, where we display the Monte Carlo sign ($), (H), (Q 2 ), (J 2 ), the 
pairing energy —g(P^P) = —g p {P^P) p — g n (P^P) n , and the moment of inertia, (X 2 ). Note 
that the sign degrades quite rapidly with increasing u, making cranking calculations at 
lower temperatures difficult. Moments of inertia were calculated from Z 2 = d(J z )/du = 
P[(J 2 ) — (Jz) 2 ]- At high temperatures, the nucleus is unpaired and the moment of inertia 
decreases as the system is cranked. However, for lower temperatures when the nucleons are 
paired, the moment of inertia initially increases as we begin to crank, but then decreases 
at larger cranking frequencies as pairs break; Figure 8.19 shows that the pairing gap also 
decreases as a function of u. It is well known that the moment of inertia depends on the 



pairing gap [93], and that initially Z 2 should increase with increasing lo. Once the pairs have 
been broken, the moment of inertia decreases. These features are evident in the figure. 

It is of particular interest to determine the quadrupole shape of a nucleus as function of 
temperature and angular momentum. It is generally expected that some nuclei may exhibit 
a sudden phase transition from a prolate to spherical shape as the temperature increases 
In addition, as the cranking frequency is increased a transition to oblate ellipsoids is 



also expected. 

In order to obtain a more detailed picture of the deformation, we examine the components 
of the quadrupole operator = r 2 Y 2 *^. Note, however, that due to rotational invariance 
of the uncranked Hamiltonian, the expectation value of any component is expected to 
vanish. On the other hand, there is an intrinsic frame for each Monte Carlo sample, in which 
it is possible to compute the three non-zero components Q' Q , Q' 2 , and Q'_ 2 (the prime is used 
to denote the intrinsic frame). The intrinsic quadrupole moments can then be related to 
the standard deformation coordinates (3 and 7 0, p8| . The task remains to compute the 
quadrupole moments in the intrinsic frame for each Monte Carlo sample. This is accom- 
plished by computing and diagonalizing the expectation value of the cartesian quadrupole 
tensor = 3xiXj — b^r 2 for each Monte Carlo sample. From the three eigenvalues, it is 
straightforward to determine the deformation parameters as in |JB| . 



Figure 8.20 shows the evolution of the shape distribution for 170 Dy at inverse temperatures 
T _1 = 0.5, 1.0, 2.0, and 3.0 MeV -1 . These contour plots show the free energy 7), 
obtained from the shape probability distribution, P(/3, 7), by 

F(/5, 7 ) = -Tln-^^ (8.21) 
p 6 sin 37 

where the j3 3 sin 37 is the metric in the usual deformation coordinates; V was obtained 
simply by binning the Monte Carlo samples in the (3 — 7 plane. As is seen from the plots, 
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deformation clearly sets in with decreasing temperature. At high temperatures, the system 
is nearly spherical, whereas at lower temperatures, especially at T -1 = 3.0 MeV -1 , there is 
a prolate minimum on the 7 = axis. 



8.9 7-soft nuclei 

Nuclei with mass number 100 < A < 140 are believed to have large shape fluctuations in 
their ground states. Associated with this softness are spectra with an approximate 0(5) 
symmetry and bands with energy spacings intermediate between rotational and vibrational. 
In the geometrical model these nuclei are described by potential energy surfaces with a 



minimum at (3 ^ but independent of 7 [97|. Some of these nuclei have been described in 



terms of a quartic five-dimensional oscillator |98|| . In the Interacting Boson Model (IBM) 
they are described by an 0(6) dynamical symmetry |99], 100|| . In the following we review 



the first fully microscopic calculations for soft nuclei with 100 < A < 140 ||101]| . 

For the two-body interaction we used a monopole (J = 0) plus quadrupole (J = 2) force 
102|1 supplemented by a collective quadrupole interaction: 



H 2 = - E v^h p LK - \x ■■ Et-rW-,, = , (8-22) 

where :: denotes normal ordering. The single particle energies and the other parameters were 
determined as described in Ref. |101 |. 



We begin discussion of our results with the probability distribution of the quadrupole mo- 
ment. The calculation of the shape distributions included the quantum- mechanical fluctua- 
tions through the variance of the Q operator for each sample, A 2 = Ti{Q 2 U a ) /Tr(U a ) — (Q) 2 . 
The shape distribution P(/3, 7) can be converted to a free energy surface as described by 
Eq. (8.19). 

The shape distributions of 128 Te and 124 Xe are shown in Fig. 8.21 at different temperatures. 
These nuclei are clearly 7-soft, with energy minima at (3 ~ 0.06 and (3 ~ 0.15, respectively. 
Energy surfaces calculated with Strutinsky-BCS using a deformed Woods-Saxon potential 
103|| also indicate 7-softness with values of (3 comparable to the SMMC values. These 



calculations predict for 124 Xe a prolate minimum with (3 ~ 0.20 which is lower than the 
spherical configuration by 1.7 MeV but is only 0.3 MeV below the oblate saddle point, and 
for 128 Te a shallow oblate minimum with (3 ~ 0.03. These 7-soft surfaces are similar to those 
obtained in the 0(6) symmetry of the IBM, or more generally when the Hamiltonian has 
mixed U(5) and 0(6) symmetries but a common 0(5) symmetry. In the Bohr Hamiltonian, 
an 0(5) symmetry occurs when the collective potential energy depends only on (3 |J7|. The 



same results are consistent with a potential energy V({3) that has a quartic anharmonicity 
98], but with a negative quadratic term that leads to a minimum at finite (3. 

The total E2 strengths were estimated from (Q 2 ) where Q = e p Q p + e n Q n is the electric 
quadrupole operator with effective charges of e p = 1.5e and e n = 0.5e, and B(E2; Of — > 2^) 
determined by assuming that most of the strength is in the 2f state. Values for B(E2; Of — > 
2f) of 663 ± 10, 2106 ± 15, and 5491 ± 36 e 2 fm 4 were found, to be compared with the exper- 
imental values |103 of 1164, 3910, and 9103 e 2 fm 4 for 124 Sn, 128 Te and 124 Xe, respectively. 



Thus, the SMMC calculations reproduce the correct qualitative trend. The 2+ excitation 
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0.0 1.0 2.0 3.0 4.0 5.0 0.0 1.0 2.0 3.0 4.0 5.0 

P P 

FIG 8.18 Canonical and grand-canonical static observables are shown for the uncranked 
170 Dy system as a function of (3. Circles represent results of canonical projection of the 
grand-canonical fields, and squares show the grand-canonical results. Lines are drawn to 
guide the eye. We show expectation values of the energy (H), the isoscalar quadrupole (Q 2 ), 

the valence nucleon number (N), and the number variance (AN) = J (N 2 ) — (N) 2 (grand- 
canonical only), the squared angular momentum (J 2 ) and the pairing energy —g(P^P), which 
is the expectation value of the pairing terms in the Hamiltonian (from |59|). 
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FIG 8.19. Grand-canonical observables for 170 Dy at various cranking frequencies and tem- 
peratures. We show the average sign ($), the isoscalar quadrupole moment (Q 2 ), the energy 
(H), the square of the angular momentum (J 2 ), the moment of inertia (X 2 ), and the expec- 
tation value of the pairing terms in the Hamiltonian, —g(P^P) (from [55|). 
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FIG 8.20. Contour plots of the free energy (as described in the text) in polar coordinates 
in the /3-7 plane are shown for the 170 Dy system. Inverse temperatures are from 0.5 (top) 
1.0, 2.0, and 3.0 MeV -1 (bottom). Contours are shown at 0.3 MeV intervals. Lighter shades 
indicate the more probable nuclear shapes (from [59f|). 
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energies were also calculated from the E2 response function. The values of 1.12 ± 0.02, 
0.96 ± 0.02 and 0.52 ± 0.01 MeV are in close agreement with the experimental values of 1.2, 
0.8 and 0.6 MeV for 124 Sn, 128 Te and 124 Xe, respectively. 

Another indication of softness is the response of the nucleus to rotations, probed by adding 
a cranking field uJ z to the Hamiltonian and examining the moment of inertia as a function 
of the cranking frequency u. For a soft nucleus one expects a behavior intermediate between 
a deformed nucleus, where the inertia is independent of the cranking frequency, and the 
harmonic oscillator, where the inertia becomes singular. This is confirmed in Fig. 8.22 which 
shows the moment of inertia X2 for 124 Xe and 128 Te as a function of u, and indicates that 
128 Te has a more harmonic character. The moment of inertia for u = in both nuclei is 
significantly lower than the rigid body value (« A3h 2 /MeV for A = 124) due to pairing 
correlations. 

Also shown in Fig. 8.22 are (Q 2 ) where Q is the mass quadrupole, the BCS-like pairing 
correlation (A' A) for the protons and ( J z ) (neutron pairing is less affected, and therefore not 
shown). Notice that the increase in I 2 as a function of u is strongly correlated with the rapid 
decrease of pairing correlations and that the peaks in I2 are associated with the onset of a 
decrease in collectivity (as measured by (Q 2 ))- This suggests band crossing along the yrast 
line associated with pair breaking and alignment of the quasi-particle spins at uj w 0.2 MeV 
((j z ) « 7h) for 128 Te and lu « 0.3 MeV ((J z ) « Ilk) for 124 Xe. The results are consistent 
with an experimental evidence of band crossing in the yrast sequence of Xe around spin 
of 10 h ||105|| . The alignment effect is clearly seen in the behavior of (J z ) at the lower 



temperature, which shows a rapid increase after an initial moderate change. Deformation 
and pairing decrease also as a function of temperature. 

The total number of J-pairs (rij = J2 a n aj) i n the various pairing channels was also 
calculated; the results for the number of correlated pairs (after subtracting the mean-field 
values) are shown in Fig. 8.23. Since the number of neutrons in 124 Xe is larger than the 
mid-shell value, they are treated as holes. For J = and J = 2 one can compare the largest 
n a j with the number of s and d bosons obtained from the 0(6) limit of the IBM. For 124 Xe 
the SMMC (IBM) results in the proton-proton pairing channel are 0.85 (1.22) s (J = 0) 
pairs, and 0.76 (0.78) d (J = 2) pairs, while in the neutron-neutron channel we find 1.76 
(3.67) s pairs and 2.14 (2.33) d pairs. For the protons the SMMC d to s pair ratio 0.89 is 
close to its 0(6) value of 0.64. However, the same ratio for the neutrons, 1.21, is intermediate 
between 0(6) and SU(3) (where its value is 1.64) and is consistent with the neutrons filling 
the middle of the shell. The total numbers of s and d pairs -1.61 proton pairs and 3.8 
neutron (hole) pairs - are below the IBM values of 2 and 6, respectively. 

8.10 Double-beta decay 

The second-order weak process (Z, A) — ► (Z + 2, A) + 2e~ + 2u e is an important "background" 
to searches for the lepton- number violating neutrinoless mode, {Z,A) — > (Z + 2, A). The 
calculation of the nuclear matrix element for these two processes is a challenging problem in 
nuclear structure, and has been done in a full pf model space for only the lightest of several 
candidates, 48 Ca. P.B. Radha et al. have performed first Monte Carlo calculations of the 



2v (5(3 matrix elements in very large model spaces [106 
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FIG 8.21. Contours of the free energy (as described in the text) in the polar-coordinate 
(3 — 7 plane for 128 Te and 124 Xe. Contours are shown at 0.3 MeV intervals, with lighter 
shades indicating the more probable nuclear shapes (from ||101|| ). 
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FIG 8.22. Observables for 124 Xe and 128 Te as a function of cranking frequency uo and for two 
temperatures. I2 is the moment of inertia, Q is the mass quadrupole moment, A is the J = 
pairing operator, and J z is the angular momentum along the cranking axis (from ||101|| ). 
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FIG 8.23. Number of correlated pairs of angular momentum J in 124 Xe. Different shades 
correspond to p-p, n-n, isoscalar and isovector p-n pairs. Neutrons are treated as holes. 



In two-neutrino double /3-decay, the nuclear matrix element of interest is 

M 2v = £ (/olGN-MGI'o) > (8 . 23) 

where |zo) and |/o) are the + ground states of the initial and final even-even nuclei, and \m) 
is a 1 + state of the intermediate odd-odd nucleus; the sum is over all such states. In this 
expression, G = crr_ is the Gamow- Teller operator for /3~-decay (i.e., that which changes 
a neutron into a proton) and Q = (Ei Q + Ef )/2. A common approximation to M 2v is the 
closure value, 

M 2u = ^£ (8.24) 
E 

where E is an average energy dominator and 

M c = 53</ |G|m)(m|G|io> = (/o|G • G|i > . (8.25) 

m 

SMMC methods can be used to calculate both M r and M 2v . To do so, consider the function 



(r,r') 



( e £(r+T') G t . GV^Ge-^'G) 



z 



(8.26) 



where Z = Tr^e _/3iT is the partition function for the initial nucleus, H is the many-body 
Hamiltonian, and the trace is over all states of the initial nucleus. The quantities (/3 — r — r') 
and t play the role of the inverse temperature in the parent and daughter nucleus respectively. 
A spectral expansion of shows that large values of these parameters guarantee cooling to 
the parent and daughter ground states. In these limits, we note that 0(r, r' = 0) approaches 
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e T ®\M C \ 2 , where Q = Ef — E® is the energy release, so that a calculation of 0(r, 0) leads 
directly to the closure matrix element. If we then define 



V (T,t)= ( T dr'^r')e-^l\ (8.27) 
Jo 

and 



M^(T,r)^^f^, (8.28) 

it is easy to see that in the limit of large r, (/?— r— r'), and T, M 2u (T, r) becomes independent 
of these parameters and is equal to the matrix element in Eq. (8.23). 
In the first applications, Radha et al. calculated the 2v matrix elements for 48 Ca and 



76 Ge ||106|| . The first nucleus allowed a benchmarking of the SMMC method against direct 
diagonalization. A large-basis shell model calculation for 76 Ge has long been waited for, as 
76 Ge is one of the few nuclei where the 2vf3f3 decay has been precisely measured and the 
best limits on the 0v decay mode have been established [|107j , |108| , |109| . 



The SMMC calculation for 48 Ca was performed for the complete pf shell using the KB3 
interaction and compared to a direct diagonalization using an implementation of the Lanc- 
zos algorithm. To avoid the sign problem (section 7), the SMMC result was obtained by 
extrapolation from a family of sign-problem-free Hamiltonian H g with g < and \ = 4. 
For both the exact and the closure matrix elements, the SMMC and direct diagonalization 
results agree well for all values of g < (Fig. 8.21). Upon extrapolation to the physical 
value g = 1 the SMMC study yields M 2v = 0.15 ± 0.07 MeV, consistent with the matrix 
element obtained by direct diagonalization (M 2l/ = 0.08 | 110 |, with the erratum in [|5T|j). 



Note that the ^-dependence was found to be linear at g < 0. However, when monitoring the 
(^-dependence for g > within the direct diagonalization, a small curvature in the extrapo- 
lation was found at g ~ 1 which could not be detected in the SMMC results at g < and so 
leads to some uncertainty in the extrapolated SMMC result. 

We note that the shell model estimate for the 48 Ca lifetime |5lj] appears to be shorter than 
the recently established experimental value ||111|| . While this experimental lifetime might 



still be compatible with complete p/-shell model calculation after slight modifications of the 
important J — 0, T — 1 and J = 1,T — matrix elements | |112| |, it appears desirable to 
calculate the 2v matrix element for 48 Ca in the complete (0dls)(lp0f) shells taking possible 
polarization effects into account. Such a challenging calculation is computationally feasible 
within the SMMC. 

To monitor the possible uncertainty related to the ^-extrapolation in the calculation of 
the 2v matrix element for 76 Ge, SMMC studies have been performed for two quite different 
families of sign-problem-free Hamiltonians (x = oo and x = 4). The calculation comprises 
the complete (O/5/2, lp, Ogg/2) model space, which is significantly larger than in previous shell 
model studies [|113j| . The adopted effective interaction is based on the Paris potential and 



has been constructed for this model space using the Q-box method developed by Kuo |114fl 



As is shown in Fig. 8.25, upon linear extrapolation both families of Hamiltonians predict a 
consistent value for the 2v matrix element of 76 Ge. The results M 2v = 0.12±0.07 and M 2v = 
0.12 ± 0.06 are only slightly lower than the experimental values [M 2v = 0.22 ± 0.01, [p9[| ). 
This comparison, however, should not be overinterpreted yet, as the detailed reliability of 
the effective interaction is still to be checked. 
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It is interesting that the closure matrix element found in the SMMC calculation and 
the average energy denominator (M c = —0.36 ± 0.37, E = —3.0 ± 3.3 MeV and M c = 
0.08 ± 0.17, E = 0.57 ± 1.26 MeV for the two families of Hamiltonians with % = oo and 
X = 4, respecitvely) are both significantly smaller than had been assumed previously. This 
is confirmed by a recent truncated diagonalization study ||115||. 



9 Prospects 

We have presented in section 8 a sampling of results from SMMC calculations. These demon- 
strate both the power and limitations of the methods and the physical insights they offer. 
SMMC calculations, while not a panacea, clearly have certain advantages over conventional 
shell model approaches, particularly for properties of ground states or thermal ensembles. 
Of the results discussed in this review, the most significant bear on the quenching of GT 
strength, the pairing structure, and nuclear shapes. Many of the first applications of the 
SMMC approach have significant bearing on astrophysical questions. 
With respect to the technical aspects of these calculations we note the following. 

• SMMC methods are computationally intensive. However, computing power is becom- 
ing both less expensive and more widely available at an astonishing rate. It is a great 
advantage that these calculations can efficiently exploit loosely connected "farms" of 
work-station-class machines. 

• A number of interesting physics questions require multi-shell model spaces. Among 
these are the origin of the apparent renormalization of qa = 1 in Gamow- Teller tran- 
sitions, strengths of first-forbidden weak operators, the properties of the giant dipole 
resonance, and nuclear behavior at temperatures above 1.5 MeV. Preliminary SMMC 
studies offer circumstantial evidence that center-of-mass (CM) motion is not a signifi- 
cant concern for some of the operators of interest. This is not too surprising at finite 
temperature, since the CM is only three degrees of freedom (far fewer than the internal 
dynamics) . 

• We lack the capability to treat odd- A or odd-odd N ^ Z systems at temperatures 
below ~800 keV, because of " sign" problems in the Monte Carlo sampling. However, 
meaningful results are possible at modest temperatures for such nuclei, as examplified 
by the 51 V, 55 Mn, and 59 Co calculations in section 8.1. Similar problems prevent spin 
projection, which would enable yrast spectroscopy. 



Otsuka and collaborators ||116|| have recently proposed a hybrid scheme whereby SMMC 
methods are used to select a many-body basis, which is then employed in a conventional 
diagonalization. The sign problems alluded to above are absent, and detailed spectroscopy 
is possible. Test applications to boson problems have shown some promise, although the 
utility for realistic fermion systems remains to be demonstrated. It is unfortunate that 
explicit angular momentum projection apparently plays an important role in these methods, 
as the numerical effort of that procedure increases strongly with the size of the model space. 

Additional physics results in finite nuclei that should emerge in the future include: more 
realistic electron capture rates in pre-supernova conditions, the double-beta decay matrix 
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FIG 8.24. The closure matrix element (upper panel) and the exact matrix element (lower 
panel) for the double-/? decay of 48 Ca as a function of g- values in the adopted family of Hamil- 
tonian (see text). The SMMC results (open circles) are compared to the direct diagonaliza- 
tion (solid circles). The SMMC result at g = 1 has been obtained by linear extrapolation 
(from MP . 
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FIG 8.25. The 2u matrix element for 76 Ge calculated within SMMC studies based on two 
families of Hamiltonians which are free of sign problems. The physical values are obtained 
by linear extrapolation to g = 1. The experimental value for this matrix element ||109|| is 
indicated by the diamond (from ||106|. 



elements for candidates other than 76 Ge, systematic studies of rare earth nuclei at finite 
temperature and spin, studies to improve the effective interactions used, tests of such models 
as the IBM and RPA, and predictions of nuclear properties far from /3-stability. 

The ground state and thermal properties of nuclear matter are another intriguing applica- 
tion of SMMC methods. One approach is to use single particle states that are plane waves 
with periodic boundary conditions and a G-matrix derived from a realistic inter-nucleon in- 
teraction; the formalism and algorithms we have presented here are then directly applicable. 
An alternative approach is to work on a regular lattice of sites in coordinate space and em- 
ploy Skyrme-like effective interactions that couple neighboring sites; the calculation is then 
similar to that for the Hubbard model for which special techniques must be used to handle 
the large, sparse matrices involved pi. Both approaches are currently being pursued. 
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A Appendix: Properties of one-fermion operators 



In this appendix, we review some of the properties of one-fermion operators relevant to Monte 



Carlo methods for the shell model. A more complete discussion can be found in Ref. |27] . 

Consider a set of N s single-particle states labeled by a = 1, . . . , N s . Corresponding to each 
state a are anticommuting fermion creation (ajj and annihilation (a a ) operators that satisfy 

{a a , ap] = {a^, ajj} = 0; {a+, a p } = 5 a/3 . (A.l) 

The associated hermitian number operators, h a = a^ a a a , are mutually commuting, and satisfy 
the operator identity 

n 2 a = n a , (A.2) 

implying that their eigenvalues are either or 1; i.e., a single-particle state can either be 
empty or occupied. The total number operator is N = }^h a . The fermion vacuum, |0), 

a 

is annihilated by all operators, a a \0) = 0, and so satisfies h a \0) = 0; i.e., all single-particle 
states are empty. 

A complete orthonormal basis of A-fermion states can be constructed by choosing A dif- 
ferent single-particle states to be occupied (let them be labeled by i = 1, . . . , A), and then 
constructing the Slater determinant 

\$) = a[a& . . . a\\0) , (A.3) 

where the product is over all occupied states and the order of the creation operators follows 
a fixed, predetermined sequence. This state is an eigenvalue of each n« with eigenvalue 1 or 
0, depending upon whether the state i is occupied or not. 
One-body operators are defined as 

O = E °^4«/3 , (A.4) 

a/3 

where the O a p = (ai\0\(3) are the one-body matrix elements. (Note that we must be careful 
to distinguish between the second-quantized operator O and the c-number matrix of its one- 
body matrix elements, O.) When acting on a determinant, a one-body operator produces 
a linear combination of other determinants that differ from the original by at most the 
occupation number of one state. 

Operators and determinants can be defined in any orthonormal single-particle basis. Sup- 
pose that there is a second basis A/i . . . distinct from the af3 . . . basis, and let T Xa = (A|a) 
be the N s x N s unitary matrix effecting transformations between them. Then the creation 
operators transform as 

4 = ]Tt Aq 4 (a.5) 

a 

and one-body matrix elements transform as 

O v = E r AA^. (A.6) 

a(3 
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Of particular relevance for SMM C are the exponentials of one-body operators. Consider the 
action of u — e~ Af3h on a determinant of the form (A. 3), where h is a one-body hamiltonian. 
Since the Baker-Hausdorff identity implies that 



e -A/V ie AM = £ [ e -A/3h] Qi at ^ X: Uai4 (A.7) 



and u\0) = |0), then 



1, 2, . . . , A) = ^£ Uaiot) u^j . . . ^ u^aj j |0) . (A.8) 



Thus, the action of the exponential of a one-body operator on a determinant is to redefine 
the occupied states by the linear transformation associated with u. For evolution at a finite 
f3, we deal with products of exponentials of one-body operators, U — . . . u^Ui] the net result 
of this operator on a determinant is clearly a linear transformation of the occupied states by 
the matrix U = . . . u 2 ui. It then follows that the expectation value of U in a determinental 
trial function, as is required in the zero-temperature formalism, is 

($\U\<S>) = destUij , (A.9) 

where Uij is the A x A matrix of one-body quantities, (i\U\j). 

The grand-canonical trace of U is the sum of the expectation values in all possible many- 
particle states: 

TtU = E ^aU , (A.10) 

A=l 

where Tr^ is the canonical trace (sum over all A-particle states). The enumeration of the 
canonical traces is straightforward: 

Tr U = (0|f)|0> = l 

TnU = EUI^I 1 ) = trU 
i 

Tr 2 t> = iE( 12 |f>|12) = ^[(trU) 2 -trU J 
z 12 z 



Tr Ns U = detU. (A.ll) 

These can be summed into the single expression 

TrU = det(l + U) , (A. 12) 

which can be verified by a direct expansion in powers of U. Here, the symbol tr denotes the 
matrix trace, while Tr is reserved for many-body traces. 
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We will also be interested in the thermodynamic expectation values of few-body operators. 
For a one-body operator Q, we require 

<fl> = ™£ . (A.13) 
TiU 

This expression is most conveniently evaluated by considering the operator U £ = Ue sn , so 
that 

<6> = ^lnTrC/ £ | £=0 . (A.14) 

Since Tr U £ = det(l + Ue e ") and det(A + eB) ~ (det A)(l + eTr A _1 B) to linear order in e, 
we have 

(Cl) = tr — !— UfJ . (A.15) 
N ' 1 + U 

Similarly the expectation value of the product of two one-body operators can be found to 
be 

(Ci'Ci) = (Ci')iCi) + tr— !— un'n-tr— ^— un'—?— un . (A.i6) 

N/NM/ 1+U 1+U 1+U v ' 
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